# Computing temperature of a fluid inside a cylinder using PDEPE

22 views (last 30 days)

Show older comments

Rodrigo Resende
on 10 Aug 2022

Commented: Bill Greene
on 15 Aug 2022

I am trying to solve the following problem:

- there is a fluid inside a cylinder with an initial temperature;
- this fluid starts cooling down while cylinder is cooled by the water surrounding it;
- 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?

##### 3 Comments

Bill Greene
on 11 Aug 2022

### Accepted Answer

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.

I have included my input file below.

function matlabAnswers_8_11_2022

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

### 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.

##### 9 Comments

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

##### 0 Comments

Rodrigo Resende
on 15 Aug 2022

Edited: Rodrigo Resende
on 15 Aug 2022

##### 1 Comment

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

add this line:

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

### See Also

### Products

### Community Treasure Hunt

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

Start Hunting!