# Temperature profile in a pipe flow with Matlab pdepe solver

24 Ansichten (letzte 30 Tage)
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.
%% 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 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Antworten (1)

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 KommentareKeine anzeigenKeine ausblenden
Julie am 15 Mai 2024
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 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.
%% 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.

### Kategorien

Mehr zu Heat Transfer finden Sie in Help Center und File Exchange

R2024a

### Community Treasure Hunt

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

Start Hunting!

Translated by