Temperature profile in a pipe flow with Matlab pdepe solver
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;
endIβ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Β 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Β fluid meganics, ohmic heating proces, heat transfer, pdepe solverΒ MATLAB Answers β New Questions
β