# Computing temperature of a fluid inside a cylinder using PDEPE

22 views (last 30 days)
Rodrigo Resende on 10 Aug 2022
Commented: Bill Greene on 15 Aug 2022
I am trying to solve the following problem:
1. there is a fluid inside a cylinder with an initial temperature;
2. this fluid starts cooling down while cylinder is cooled by the water surrounding it;
3. my goal is to plot "time vs. fluid temperature".
I started modeling the problem using "pdepe" funcion. What I did first was trying to model just the region between R0=4in and R1=6in and trying to apply the boundary conditions to the inner and outer faces. I had no problems with the outer face, because water outside cylinder has a fixed temperature and convection coefficient is also a constant.
However, I didn't manage to model the inner boundary condition. For the purposes of what I need, all the fluid inside the cylinder can be consider to have the same temperature at a given time. Therefore, it can be considered that , where T is the the fluid temperature inside the cylinder and is the temperature of the inner wall of the cylinder (it is a funcion of the time and will be computed by PDEPE).
I did not manage to write this BC as it is required by PDEPE function, .
Anyone can help?
Bill Greene on 11 Aug 2022
Yes,I was thinking of trying pde1dm on this problem (using one ODE as you suggest) if I could understand the physics. ;-)

Bill Greene on 12 Aug 2022
I have written a 1D PDE solver that has an input syntax similar to pdepe but has some
additional enhancements including the capability to couple additional ODEs to the
PDE. (This is the code Torsten mentioned in his comments.)
I have used this code to model your problem, as I understand it. Since you didn't provide
any values for the material and other properties, I have simply assumed some that seem
reasonable to me. The computed results also look reasonable to me. The code has been
well-tested so I'm fairly confident of its correctness but I am less confident I have
correctly modeled your particular problem.
If you want to try this example yourself, the solver can be downloaded here.
I have included my input file below.
R0=.1; R1=.15; % radii in meters
wallConductivity=150;
wallDensity=2800;
wallSpecHeat=800;
hWater=100; % convection coeff between tank and containing water
hTank=200; % convection coeff between tank and inner reservoir
tempTank=50; % initial temperature of fluid in the inner reservoir
tempWater=20; % temperature of containing water
x = linspace(R0,R1,10);
t = linspace(0,1e4,50);
pdeFunc = @(x,t,u,DuDx) thepde(x,t,u,DuDx,wallConductivity,wallDensity,wallSpecHeat);
initialWallTemp = tempWater;
icFunc = @(x) theIC(x, initialWallTemp);
bcFunc = @(xl,ul,xr,ur,t, v, vdot) theBC(xl,ul,xr,ur,t,v, vdot, hWater,tempWater,hTank);
odef = @(t,v,vdot,x,u) odeFunc(t,v,vdot,x,u,hTank);
odeicf=@() odeIcFunc(tempTank);
opts.abstol=1e-8;
opts.reltol=1e-4;
m=1;
xOde=R0;
[u,uOde] = pde1dm(m, pdeFunc,icFunc,bcFunc,x,t,odef, odeicf,xOde, opts);
figure; plot(t, u(:,1)); grid; xlabel t; title('Temperature at inner surface');
figure; plot(t, u(:,end)); grid; xlabel t; title('Temperature at outer surface');
figure; plot(t, uOde); grid; xlabel t; title 'Tank temperature';
end
function [c,f,s] = thepde(x,t,u,DuDx,wallConductivity,wallDensity,wallSpecHeat)
c = wallDensity*wallSpecHeat;
f = wallConductivity*DuDx;
s=0;
end
function u0 = theIC(x, initialWallTemp)
u0 = initialWallTemp;
end
function [pl,ql,pr,qr] = theBC(xl,ul,xr,ur,t,v, vdot, hWater, tempWater, hTank)
pl = hTank*(v-ul);
ql = 1;
pr = hWater*(ur-tempWater);
qr = 1;
end
function f=odeFunc(t,v,vdot,x,u,hTank)
% the ODE variable defines the temperature of the fluid in the tank
rho=997;
cp=4182;
R0=x;
f=vdot+2*hTank*(v-u)/(rho*cp*R0);
end
function v0=odeIcFunc(tempTank)
v0=tempTank;
end
##### 1 CommentShowHide None
Rodrigo Resende on 12 Aug 2022
Thank you, @Torsten and @Bill Greene!
I'll try to use your solver this weekend (or today) and then I tell you what I get.

### More Answers (3)

Torsten on 10 Aug 2022
Therefore, it can be considered that , Tdot = -2*h*(T-Tw)/(rho*cp*R0) where T is the the fluid temperature inside the cylinder and Tw is the temperature of the inner wall of the cylinder (it is a funcion of the time and will be computed by PDEPE).
And how do you get T ? Or is this your question ?
If you have T available, the equation for pdepe at r = r0 would be
lambda*dTw/dr = alpha*(Tw-T)
This can be set in bcfun as
pl = -alpha*(ul-T);
ql = 1.0;
if you have set
f = lambda*dudx
in pdefun.
Torsten on 12 Aug 2022
Edited: Torsten on 12 Aug 2022
I think you can get a reliable solution if you call "pdepe" in a loop from t0 to t1, from t1 to t2 and so on. During the time intervals t_i <= t <= t_i+1, you can keep T constant at the value you received with your Euler forward at time t_i.
Or try Bill Greene's extended version of pdepe (pde1dm):
I think the code can handle additional ODEs.
Note that the additional x under the operator
1/x * d/dx (x*du/dx)
is not specified when you define the boundary condition for pdepe.
Thus your pl,ql,pr,qr in "heatbc" are wrong.

Torsten on 12 Aug 2022
Edited: Torsten on 12 Aug 2022
For comparison:
R0 = .1;
R1 = .15; % radii in meters
t0 = 0;
t1 = 1e4; % time limits
wallConductivity = 150;
wallDensity = 2800;
wallSpecHeat = 800;
rhoWater = 997;
cpWater = 4182;
hWater = 100; % convection coeff between tank and containing water
hTank = 200; % convection coeff between tank and inner reservoir
tempTank = 50; % initial temperature of fluid in the inner reservoir
tempWater = 20; % temperature of containing water
nr = 10;
r = linspace(R0,R1,nr); %spatial discretization
r = r.';
rm = 0.5*(r(1:nr-1)+r(2:nr)); %cell centers
nt = 50;
t = linspace(t0,t1,nt); %output times
initialtempWall = tempWater; %initial temperatures
initialtempTank = tempTank;
y0 = [initialtempWall*ones(nr,1);initialtempTank];
[T,Y] = ode15s(@(t,y)odefun(t,y,nr,r,rm,wallConductivity,wallDensity,wallSpecHeat,...
rhoWater,cpWater,hWater,hTank,tempWater),t,y0); % solver call
%Plot temperature profile in the wall for certain output times
figure(1)
plot(r,[Y(1,1:nr);Y(2,1:nr);Y(5,1:nr);Y(10,1:nr);Y(11,1:nr)]) %Plot temperature profiles in certain wall points over time
figure(2)
plot(t,[Y(:,1),Y(:,5),Y(:,10)]) %Plot temperature of water in inner cylinder
figure(3)
plot(t,Y(:,nr+1)) function dTemp = odefun(t,y,nr,r,rm,wallConductivity,wallDensity,wallSpecHeat,...
rhoWater,cpWater,hWater,hTank,tempWater)
tempWall = y(1:nr);
tempTank = y(nr+1);
dtempWalldt = zeros(nr,1);
% Wall Temperatures
% left boundary
dtempWalldt(1) = 1/r(1)*(rm(1)*wallConductivity*(tempWall(2)-tempWall(1))/(r(2)-r(1))-...
r(1)*hTank*(tempWall(1)-tempTank))/(rm(1)-r(1))/...
(wallDensity*wallSpecHeat);
% inner points
dtempWalldt(2:nr-1) = 1./r(2:nr-1).*(rm(2:nr-1).*wallConductivity.*(tempWall(3:nr)-tempWall(2:nr-1))./(r(3:nr)-r(2:nr-1)) -...
rm(1:nr-2).*wallConductivity.*(tempWall(2:nr-1)-tempWall(1:nr-2))./(r(2:nr-1)-r(1:nr-2)))./(rm(2:nr-1)-rm(1:nr-2))./...
(wallDensity*wallSpecHeat);
% right boundary
dtempWalldt(nr) = 1/r(nr)*(r(nr)*hWater*(tempWater-tempWall(nr))-rm(nr-1)*wallConductivity*...
(tempWall(nr)-tempWall(nr-1))/(r(nr)-r(nr-1)))/(r(nr)-rm(nr-1))/...
(wallDensity*wallSpecHeat);
% Tank Temperature
dtempTankdt = 2/r(1) * hTank* (tempWall(1) - tempTank)/(rhoWater*cpWater);
% Return temporal derivatives
dTemp = [dtempWalldt;dtempTankdt];
end

Rodrigo Resende on 15 Aug 2022
Edited: Rodrigo Resende on 15 Aug 2022
I have addapted the script @Bill Greene had sent and it seems perfect to me. As long as in my problem the wall of the cylinder has many layers with different materials, I included some "if statements" to decide which properties to use as a function o "x" (did that inside "thepde" function that Bill sent).
I have compared the results with the ones obtained with a comercial simulator and it seems just the same. A reference fluid cooling from T1=55°C down to T2=25°C took circa 6.8h in both software, with just 1.25min gap between them. And it may be even closer, because I simplified the internal convective heat coefficient, the density and the specific heat of the internal fluid in order to run the matlab model.
So, I'm pretty happy with the results and the flexibility that your script has. It seems a really powerful tool for these type of problems. Thank you, guys!
As I ran in an old version of Matlab (2010b), I had to do the following modifications in the code I got from GitHub (I put it here because it may help others):
2. PDE1dImpl.m: "At least one of the entries in the c-coefficient vector must be non-zero." ---> 'At least one of the entries in the c-coefficient vector must be non-zero.'
3. PDE1dImpl.m: 'like', depVarClassType ---> class(depVarClassType)
dNdx = dN(:,i)./jac ----> dNdx = dN(:,i)./jac(1)
4. PDEMeshMapper.m: 'like', srcU ----> class(srcU)
The only change I was in doubt was if I could replace "jac" by "jac(1)".
##### 1 CommentShowHide None
Bill Greene on 15 Aug 2022
Glad to hear that pde1dm worked for you!
>The only change I was in doubt was if I could replace "jac" by "jac(1)".
Apparently the semantics of the ./ operator were enhanced sometime after the 2010b version.
The change you made is OK as long as all the elements in your mesh (the entries in the xmesh array)
are uniformly spaced. But, for a non-uniform mesh, it will fail because the values
in the jac array vary.
A change that should work in general, possibly at the cost of being slower, is to
dNdx=zeros(nen, length(jac));
after this line:
jac=(x(i2)-x(i1))/2;
and replace
dNdx = dN(:,i)./jac;
with
for j=1:nen
dNdx(j,:) = dN(j,i)./jac;
end