清空筛选条件
清空筛选条件

Temperature profile in a pipe flow with Matlab pdepe solver

6 次查看（过去 30 天）
Julie
Julie 2024-5-15
编辑： Torsten 2024-5-16
Ran in:
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

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2024-5-15
编辑：Torsten 2024-5-15
Ran in:
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 个评论
Julie
Julie 2024-5-15
thankyou for replying and adjusting!
How do you mean my Q is different from the article?
Q= elec_cond * |▽V|^2
=> elec_cond * E^2
with E= V /L = (voltage)/ (length of the pipe)
Torsten
Torsten 2024-5-15
编辑：Torsten 2024-5-16
Ran in:
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Heat Transfer 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by