CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Convection-diffusion energy equation with specified heat transfer

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree10Likes
  • 1 Post By arjun
  • 1 Post By sbaffini
  • 1 Post By mcanoenen
  • 1 Post By sbaffini
  • 1 Post By sbaffini
  • 1 Post By sbaffini
  • 1 Post By sbaffini
  • 1 Post By FMDenaro
  • 1 Post By sbaffini
  • 1 Post By sbaffini

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 14, 2021, 21:51
Default Convection-diffusion energy equation with specified heat transfer
  #1
New Member
 
Umberto
Join Date: Sep 2021
Posts: 6
Rep Power: 2
UmbertoBitencourt is on a distinguished road
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.

Code:
%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.
Attached Images
File Type: jpg Comparison.jpg (28.5 KB, 11 views)
UmbertoBitencourt is offline   Reply With Quote

Old   September 15, 2021, 00:44
Default
  #2
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,079
Rep Power: 29
arjun will become famous soon enougharjun will become famous soon enough
Try smaller dt than what you are using, it may happen that you time step size is much larger than what shall be.
UmbertoBitencourt likes this.
arjun is offline   Reply With Quote

Old   September 15, 2021, 04:33
Default
  #3
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 1,847
Blog Entries: 29
Rep Power: 36
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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.
UmbertoBitencourt likes this.
sbaffini is offline   Reply With Quote

Old   September 15, 2021, 04:48
Default Error limit
  #4
New Member
 
Murat Can
Join Date: Sep 2021
Location: Turkey
Posts: 8
Rep Power: 2
mcanoenen is on a distinguished road
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
UmbertoBitencourt likes this.
mcanoenen is offline   Reply With Quote

Old   September 15, 2021, 05:10
Default
  #5
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 1,847
Blog Entries: 29
Rep Power: 36
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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 likes this.
sbaffini is offline   Reply With Quote

Old   September 15, 2021, 08:45
Default
  #6
New Member
 
Umberto
Join Date: Sep 2021
Posts: 6
Rep Power: 2
UmbertoBitencourt is on a distinguished road
Quote:
Originally Posted by sbaffini View Post
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 is offline   Reply With Quote

Old   September 15, 2021, 08:49
Default
  #7
New Member
 
Umberto
Join Date: Sep 2021
Posts: 6
Rep Power: 2
UmbertoBitencourt is on a distinguished road
Quote:
Originally Posted by mcanoenen View Post
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.
UmbertoBitencourt is offline   Reply With Quote

Old   September 15, 2021, 08:49
Default
  #8
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 1,847
Blog Entries: 29
Rep Power: 36
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Quote:
Originally Posted by UmbertoBitencourt View Post
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?
UmbertoBitencourt likes this.
sbaffini is offline   Reply With Quote

Old   September 15, 2021, 09:00
Default
  #9
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 1,847
Blog Entries: 29
Rep Power: 36
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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 likes this.
sbaffini is offline   Reply With Quote

Old   September 15, 2021, 09:53
Default
  #10
New Member
 
Umberto
Join Date: Sep 2021
Posts: 6
Rep Power: 2
UmbertoBitencourt is on a distinguished road
Quote:
Originally Posted by sbaffini View Post
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.
UmbertoBitencourt is offline   Reply With Quote

Old   September 15, 2021, 10:02
Default
  #11
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 1,847
Blog Entries: 29
Rep Power: 36
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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 likes this.
sbaffini is offline   Reply With Quote

Old   September 15, 2021, 10:15
Default
  #12
New Member
 
Umberto
Join Date: Sep 2021
Posts: 6
Rep Power: 2
UmbertoBitencourt is on a distinguished road
Quote:
Originally Posted by sbaffini View Post
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.
UmbertoBitencourt is offline   Reply With Quote

Old   September 15, 2021, 12:06
Default
  #13
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,099
Rep Power: 66
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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.
UmbertoBitencourt likes this.
FMDenaro is online now   Reply With Quote

Old   September 15, 2021, 12:12
Default
  #14
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 1,847
Blog Entries: 29
Rep Power: 36
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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):

\frac{\partial \left(UT\right)}{\partial x} = \frac{\partial}{\partial x}\left[\left(\alpha + \alpha_t \right)\frac{\partial T}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\alpha + \alpha_t \right)\frac{\partial T}{\partial y}\right]

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 \alpha_t, 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 likes this.
sbaffini is offline   Reply With Quote

Old   September 16, 2021, 00:31
Default
  #15
New Member
 
Umberto
Join Date: Sep 2021
Posts: 6
Rep Power: 2
UmbertoBitencourt is on a distinguished road
Quote:
Originally Posted by sbaffini View Post
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):

\frac{\partial \left(UT\right)}{\partial x} = \frac{\partial}{\partial x}\left[\left(\alpha + \alpha_t \right)\frac{\partial T}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\alpha + \alpha_t \right)\frac{\partial T}{\partial y}\right]

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 \alpha_t, 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.
UmbertoBitencourt is offline   Reply With Quote

Old   September 16, 2021, 03:18
Default
  #16
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 1,847
Blog Entries: 29
Rep Power: 36
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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
UmbertoBitencourt likes this.
sbaffini is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Coupled Heat and Mass Transfer Mecroob OpenFOAM Running, Solving & CFD 1 July 12, 2020 19:24
oscillation of energy equation in steady state heat transfer mohammadreza_hj Fluent UDF and Scheme Programming 0 June 30, 2020 11:36
Interphase mass transfer of a reaction cfx_ws1992 Main CFD Forum 0 May 15, 2017 21:42
ATTENTION! Reliability problems in CFX 5.7 Joseph CFX 14 April 20, 2010 15:45
Two-Phase Buoyant Flow Issue Miguel Baritto CFX 4 August 31, 2006 12:02


All times are GMT -4. The time now is 14:39.