[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.

[arjun]

Try smaller dt than what you are using, it may happen that you time step size is much larger than what shall be.

[sbaffini]

What scheme are you using?

Also, how did you get the "ANSYS" reference solution? I'm pretty sure it's not straightforward at all, given the velocity profile you are using.

Finally, with prescribed velocity you could just solve a single linear system once, instead of iterating.

[mcanoenen]

Hello, have you tried to take the error limit like err=10^-10, sometimes, even though the error is so small you can encounter an nonconvergent solution when you check a property such as temperature in your case. In my study I am working on my CFD code (which is very simplific). I check the residuals (sum of the errors in each cell), and I can see even 10^-4 residual can cause my flow field change as the number of iterations incerases.

Regards

[sbaffini]

Not directly related to the issue but, if you don't want to go implicit, at least precompute the whole alfa, no need to do that in a loop.

[UmbertoBitencourt]

Originally Posted by sbaffini

What scheme are you using?

Also, how did you get the "ANSYS" reference solution? I'm pretty sure it's not straightforward at all, given the velocity profile you are using.

Finally, with prescribed velocity you could just solve a single linear system once, instead of iterating.

I solved the problem on ANSYS before and exported the solution to MATLAB.
Yeah, you are right, I can solve for the alfa's outside of the loop. Thanks man!

[UmbertoBitencourt]

Originally Posted by mcanoenen

Hello, have you tried to take the error limit like err=10^-10, sometimes, even though the error is so small you can encounter an nonconvergent solution when you check a property such as temperature in your case. In my study I am working on my CFD code (which is very simplific). I check the residuals (sum of the errors in each cell), and I can see even 10^-4 residual can cause my flow field change as the number of iterations incerases.

Regards

I tried that too. It didn't work. Thanks for your answer anyway.

[sbaffini]

Originally Posted by UmbertoBitencourt

I solved the problem on ANSYS before and exported the solution to MATLAB.
Yeah, you are right, I can solve for the alfa's outside of the loop. Thanks man!

What I meant is that you are using a very specific analytical velocity profile. Are you sure that the solution you got in ANSYS (Fluent?) is compatible with this velocity profile in the first place? For example, this would have required you to first initialize the velocity field with your analytical solution and then solve the temperature without solving the flow, assuming all the mass fluxes etc. were correctly precomputed. Is this what you did?

[sbaffini]

Also, I'm pretty sure that there is no turbulence model similar to the VanDriest one in any ANSYS CFD product (yet, SA might sometime resemble it). So, you are definitely not solving the same problem as in ANSYS.

So, before attempting this, I would definitely first try an equivalent laminar problem, with parabolic velocity profile and no alpha_t

[UmbertoBitencourt]

Originally Posted by sbaffini

Also, I'm pretty sure that there is no turbulence model similar to the VanDriest one in any ANSYS CFD product (yet, SA might sometime resemble it). So, you are definitely not solving the same problem as in ANSYS.

So, before attempting this, I would definitely first try an equivalent laminar problem, with parabolic velocity profile and no alpha_t

Exactly, there is no equivalent on ANSYS. But I expected that the answer would make sense at least. I even tried a constant velocity profile with no alfa_t and the answer was pretty much the same: that flat plot that looks like a step function.

I was wondering if I should have included the alfa*dT^2/dx^2 in the PDE. But that would require one more BC and I don't think I can provide one more.

[sbaffini]

Without knowing, exactly, which scheme you are using and how the case was setup in Fluent, it's impossible to help.

For example, I now seem to understand that you used the x direction as a sort of (pseudo)time?

If this is the case, and you're not really interested in the x evolution, you are doing this all wrong. You would just need a single TDMA solve (1D).

So, may I ask what you want to achieve and from what input (what started you doing this in the first place)?

[UmbertoBitencourt]

Originally Posted by sbaffini

Without knowing, exactly, which scheme you are using and how the case was setup in Fluent, it's impossible to help.

For example, I now seem to understand that you used the x direction as a sort of (pseudo)time?

If this is the case, and you're not really interested in the x evolution, you are doing this all wrong. You would just need a single TDMA solve (1D).

So, may I ask what you want to achieve and from what input (what started you doing this in the first place)?

Okay, I got you. Let me try to explain better.

I'm looking for the temperature profile between two parallel plates T(x,y). This is basically simulating the test section of a 2D wind tunnel which is being heated, where I know what are the conditions at the inlet (velocity and temperature).

Do you think I might be solving the wrong eq.? I'm not interested in the transient, only in the steady-state.

I really appreciate your time, thanks for trying to understand.

[FMDenaro]

I don't understand your case, if T=T(y) the the LHS of your equation is simply zero and you have a standard tridiagonal system to solve in an exact way.

[sbaffini]

Ok, your problem seems to be the 2D channel flow, and you are interested in the stramwise (x) as well as wall normal (y) description, so you have an inlet, an outlet and two walls.

You have an assigned velocity profile U, which is constant in x. Then, the equation you should be trying to solve is (assuming constant density and Cp):



So, yes, the first thing is to include the streamwise diffusion, there is no reason to not have that but keep the convection along x when U does not depend on x. Should be evident from the equation.

Second thing, how are the terms discretized? Finite difference? Finite volume? This may have implications on how to set the bcs. However, very generally, you just need T at inlet (you already know U and can compute dT/dx from interior) and set outflow like bc at outlet (extrapolate T from interior, dT/dx = 0). At walls, you either know T and can compute dT/dy from interior, or viceversa.

You need to be careful in how you discretize the y diffusion term for variable , but the 1st order upwind is ok for the convection term for now.

Finally, once you discretize it, you get a system of equations. You can use SOR on this system (as I guess you were already doing) by expressing the T(i,j) unknown as function of the other temperatures around.

Is this what you're doing?

[UmbertoBitencourt]

Originally Posted by sbaffini

Ok, your problem seems to be the 2D channel flow, and you are interested in the stramwise (x) as well as wall normal (y) description, so you have an inlet, an outlet and two walls.

You have an assigned velocity profile U, which is constant in x. Then, the equation you should be trying to solve is (assuming constant density and Cp):



So, yes, the first thing is to include the streamwise diffusion, there is no reason to not have that but keep the convection along x when U does not depend on x. Should be evident from the equation.

Second thing, how are the terms discretized? Finite difference? Finite volume? This may have implications on how to set the bcs. However, very generally, you just need T at inlet (you already know U and can compute dT/dx from interior) and set outflow like bc at outlet (extrapolate T from interior, dT/dx = 0). At walls, you either know T and can compute dT/dy from interior, or viceversa.

You need to be careful in how you discretize the y diffusion term for variable , but the 1st order upwind is ok for the convection term for now.

Finally, once you discretize it, you get a system of equations. You can use SOR on this system (as I guess you were already doing) by expressing the T(i,j) unknown as function of the other temperatures around.

Is this what you're doing?

Thanks for your effort in helping me out. I took some time to reply because I was trying the approach is suggested. The results I'm obtaining right now make much more sense, but some tuning is necessary.

Answering your questions: I'm using finite volume method. I'm doing exactly as you said.

By using your approach another boundary condition is needed (as you mentioned). I tried saying that dT/dx=0 at x=L but it turned out that the maximum temperature was higher than expected.

By using the conservation of energy, the average temperature at the exit must be:
T(L, AVERAGE)=T_0+q_0*L/mC_p. Where m is the mass flow rate.

I'm getting an average temperature much higher than that at the outlet. I guess it's because of the BC at the exit.

You said I could also extrapolate it. How could I do that??

Thanks a lot for your help.

[sbaffini]

Well, as you actually have velocity and temperature profiles, the definition is not exactly that one for the average temlerature at outlet. You have that the difference between the mass averaged temperatures is q/Cp, but the velocity has to be included in the integrals of T.

Besides this, a FV approach simply requires:

1) U and T at inlet, that you know from the profiles. These enter directly in the convective term and T can be used for the diffusion term with dT/dx = [T(1,j)-Tinlet]/(dx/2)

2) Set 0 (don't set it at all) for the diffusion term at outlet. This is the extrapolation (to 1st order), assuming that Toutlet(j) = T(nx,j) so that the strramwise diffusion term turns out to be 0