Temperature profile in a pipe flow with Matlab pdepe solver
2 次查看(过去 30 天)
显示 更早的评论
I’m currently trying to solve the heat equation for the Ohmic heating of a pipe flow using the MATLAB pdepe solver. The problem is described in the following paper: https://doi.org/10.4236/ojmsi.2021.91002. This paper also contains the solution for the axial bulk temperature profile. Unfortunately, I don't find the same solution, and I wonder if anyone could help me find the error.
The simplified heat equation for this case becomes:
𝜌𝐶𝑝𝑣𝑧∂𝑇∂𝑧=𝑘𝑟∂∂𝑟(𝑟∂𝑇∂𝑟)+𝑄̂.𝑄̂=𝜎𝐸2 is the heat source due to the Ohmic heating (this is also defined in the paper). 𝑣𝑧=2𝑉𝑎𝑣𝑔(1−(𝑥𝐷/2)2) is the classical solution for the axial velocity of the fluid in the pipe.
Please find the values of the constants in the code below.
The boundary conditions are:
𝑇(𝑧=0,𝑟)=40°𝐶,∂𝑇∂𝑟|𝑧,𝑟=0=0,−𝑘∂𝑇∂𝑟|𝑧,𝑟=𝐷/2=ℎ(𝑇(𝑧,𝑟=𝐷/2)−𝑇∞).
Please find my Matlab code below. Note that in order to be able to use the Matlab pdepe solver, the time argument in the solver corresponds with the axial coordinate 𝑧. Another important remark is that I want to solve it with the option 𝑚=0.
Any help/comments/remarks will be greatly appreciated. Thanks in advance!
%% Main code
m=0;
r=linspace(0,0.025,1000); % m r corresponds with x
z=linspace(0,2.5,1000); % m z corresponds with t
sol=pdepe(m,@pdex1pde,@ic1pde,@bc1pde,r,z);
u=sol(:,:,1); % u=T
% Figure I'm trying to recreate
figure
plot(z,u(:,1))
title('Axial bulk temperature distribution')
xlabel('z')
ylabel('T(z,r=0)')
%% Functions
% Definition of the PDE
function [c,f,s]=pdex1pde(x,t,u,DuDx)
% Parameters
T_0=40; %°C
k=0.6178; % W/m*°C
rho=992.2; % kg/m^3
C_p=4182; % J/kg*°C
m=82.5/3600; % kg/s
R=25*10^(-3); % m
D=50*10^(-3); % m
A=pi*R^2; % m^2
V_gem=m/(rho*A); % kg m
k_0=0.015; % 1/°C^3 /s
elec_cond=0.998*(1+k_0*(u-T_0)); % S/m
v_z=2*V_gem*(1-(x^2/(D/2)^2)); % kg m^3 /s
Volt=2304; % V
L=2.5; % m
Q=elec_cond*(Volt/L)^2*x/k;
c=(rho*C_p*v_z*x/k);
f=x*DuDx;
s=Q;
end
% "Initial" condition
function u0=ic1pde(x)
u0=40;
end
% Boundary conditions
function [pl,ql,pr,qr]= bc1pde(xl,ul,xr,ur,t)
% Parameters
T_amb=20; % °C
T_0=40; % °C
h=10; % W/m^2*°C
k=0.6178; % W/m*°C
pl=0; % left is ul=0
ql=1; % p+qf=0
pr=h*(ur-T_amb);
qr=k/xr;
end
0 个评论
回答(1 个)
Torsten
2024-5-15
编辑:Torsten
2024-5-15
I corrected as much as I could, but the computation of Q looks different in your code in comparison to the article.
%% Main code
m=1;
r=linspace(0,0.025,1000); % m r corresponds with x
z=linspace(0,2.5,1000); % m z corresponds with t
sol=pdepe(m,@pdex1pde,@ic1pde,@bc1pde,r,z);
u=sol(:,:,1); % u=T
% Figure I'm trying to recreate
figure
plot(z,u(:,1))
title('Axial bulk temperature distribution')
xlabel('z')
ylabel('T(z,r=0)')
%% Functions
% Definition of the PDE
function [c,f,s]=pdex1pde(x,t,u,DuDx)
% Parameters
T_0=40; %°C
k=0.6178; % W/m*°C
rho=992.2; % kg/m^3
C_p=4182; % J/kg*°C
m=82.5/3600; % kg/s
R=25*10^(-3); % m
D=50*10^(-3); % m
A=pi*R^2; % m^2
V_gem = m/rho/(2*pi*integral(@(r)(1-(r/R).^2).*r,0,R));
v_z = V_gem*(1-(x/R)^2);
%V_gem=m/(rho*A); % kg m
k_0=0.015; % 1/°C^3 /s
elec_cond=0.998*(1+k_0*(u-T_0)); % S/m
%v_z=2*V_gem*(1-(x^2/(D/2)^2)); % kg m^3 /s
Volt=2304; % V
L=2.5; % m
Q=elec_cond*(Volt/L)^2;
c=rho*C_p*v_z;
f=k*DuDx;
s=Q;
end
% "Initial" condition
function u0=ic1pde(x)
u0=40;
end
% Boundary conditions
function [pl,ql,pr,qr]= bc1pde(xl,ul,xr,ur,t)
% Parameters
T_amb=20; % °C
T_0=40; % °C
h=10; % W/m^2*°C
pl=0; % left is ul=0
ql=1; % p+qf=0
pr=h*(ur-T_amb);
qr=1;
end
2 个评论
Torsten
2024-5-15
编辑:Torsten
2024-5-16
Equation (9) of the article looks strange.
The left-hand side should be rho*cp*v_z(r)*dT/dz, not rho*cp*2*W*(1-(2*r/D)^2)*dT/dz with W being the mass flow rate. The units don't work out correctly.
Further, the density of water in Table 3 is given as 992 g/m^3 - also wrong.
And don't think the temperature profile is flat over the radius such that you only need to plot T at r = 0. You must compute a kind of mean temperature over the cross sections at the various z-positions like
sol=pdepe(m,@pdex1pde,@ic1pde,@bc1pde,r,z);
u=sol(:,:,1); % u=T
% Compute mean temperature over cross sections
for i = 1:numel(z)
Tavg(i) = 2/r(end)^2*trapz(r,u(i,:).*r);
end
Some people even use mass-weighted averaging of temperature by taking the actual velocity at the r-position into account.
This doesn't look so bad:
%% Main code
m=1;
r=linspace(0,0.025,1000); % m r corresponds with x
z=linspace(0,2.5,1000); % m z corresponds with t
sol=pdepe(m,@pdex1pde,@ic1pde,@bc1pde,r,z);
u=sol(:,:,1); % u=T
% Compute mean temperature over cross sections (Area-weighted)
for i = 1:numel(z)
Tavg_area(i) = 2/r(end)^2*trapz(r,u(i,:).*r);
end
% Compute mean temperature over cross sections (Mass-weighted)
rho=992.2; % kg/m^3
m=82.5/3600; % kg/s
R=25*10^(-3); % m
V_gem = m/rho/(2*pi*integral(@(r)(1-(r/R).^2).*r,0,R));
v_z = @(r)V_gem*(1-(r/R).^2);
den = integral(@(r)v_z(r).*r,0,R);
for i = 1:numel(z)
Tavg_mass(i) = trapz(r,v_z(r).*u(i,:).*r);
end
Tavg_mass = Tavg_mass/den;
hold on
plot(z,Tavg_area,'b')
plot(z,Tavg_mass,'r')
hold off
grid on
%% Functions
% Definition of the PDE
function [c,f,s]=pdex1pde(x,t,u,DuDx)
% Parameters
T_0=40; %°C
k=0.6178; % W/m*°C
rho=992.2; % kg/m^3
C_p=4182; % J/kg*°C
m=82.5/3600; % kg/s
R=25*10^(-3); % m
D=50*10^(-3); % m
A=pi*R^2; % m^2
V_gem = m/rho/(2*pi*integral(@(r)(1-(r/R).^2).*r,0,R));
v_z = V_gem*(1-(x/R)^2);
%V_gem=m/(rho*A); % kg m
k_0=0.015; % 1/°C^3 /s
elec_cond=0.998*(1+k_0*(u-T_0)); % S/m
%v_z=2*V_gem*(1-(x^2/(D/2)^2)); % kg m^3 /s
Volt=2304; % V
L=2.5; % m
Q=elec_cond*(Volt/L)^2;
c=rho*C_p*v_z;
f=k*DuDx;
s=Q;
end
% "Initial" condition
function u0=ic1pde(x)
u0=40;
end
% Boundary conditions
function [pl,ql,pr,qr]= bc1pde(xl,ul,xr,ur,t)
% Parameters
T_amb=20; % °C
T_0=40; % °C
h=10; % W/m^2*°C
pl=0; % left is ul=0
ql=1; % p+qf=0
pr=h*(ur-T_amb);
qr=1;
end
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!