Filter löschen
Filter löschen

Temperature profile in a pipe flow with Matlab pdepe solver

24 Ansichten (letzte 30 Tage)
Julie
Julie am 15 Mai 2024
Bearbeitet: Torsten am 16 Mai 2024
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

Antworten (1)

Torsten
Torsten am 15 Mai 2024
Bearbeitet: Torsten am 15 Mai 2024
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 Kommentare
Julie
Julie am 15 Mai 2024
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 am 15 Mai 2024
Bearbeitet: Torsten am 16 Mai 2024
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

Melden Sie sich an, um zu kommentieren.

Produkte


Version

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by