Computing temperature of a fluid inside a cylinder using PDEPE

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]

I find this question very confusing. My impression is that you have some fluid in a cylinder, from r=0 up to r=R0, that is surrounded by a solid wall from r=R0 to r=R1. This solid wall is surrounded by a large volume of water that can be assumed to have a fixed temperature. But it appears you are assuming that the fluid inside ( r=0 up to r=R0) is being cooled by the containing solid cylinder but for some reason still has a uniform temperature throughout? How can this be?

[Torsten]

@Bill Greene

If R0 is small and maybe the fluid inside is stirred, temperature assimilates quite fast.

But by the way: isn't your pdepe code capable to solve this problem (one coupled ODE at the left boundary point) ?

[Bill Greene]

Yes,I was thinking of trying pde1dm on this problem (using one ODE as you suggest) if I could understand the physics. ;-)

Answer 1

[Bill Greene]

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.

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

[Rodrigo Resende]

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.

Answer 2

[Torsten]

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.

[Rodrigo Resende]

Thank you Torsten, but that's the point: I do not have T. T here is the temperature of the fluid inside the cylinder and will change its temperature as it exchange heat with the inner wall of the cylinder.

The complete system is a cylinder filled with a fluid (from r=0 up to r=R0), then the wall of this cylinder (from r=R0 up to r=R1). The whole cylinder is immersed in water. I want to find the behavior of the fluid temperature inside the cylinder, which I called T.

[Torsten]

You could define small time steps in tspan for pdepe and call an OutputFcn function in which you update T according to the wall temperature Tw from pdepe (e.g. with simple Euler forward). Then for the next time step, you could use this T to define the boundary condition at the inner wall for pdepe.

But I think a simple discretization of the cylindrical wall and a solution of the usual

rho_w*cp_w*dT_w/dt = 1/r d/dr(r*lambda_w*dT_w/dr)

using the method-of-lines (without pdepe) should be more flexible and simple.

And here it is easy to include an additional ODE for T.

[Rodrigo Resende]

I think it will work. I'll try it later today and give you a feedback.

[Rodrigo Resende]

@Torsten, could you please show me a way of using the OutputFcn to do this job? I'm not used to it and I really do not understand how it works.

I tried to do like this:

options = odeset('OutputFcn',@Tfluid);

sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,tspan,options);

As I read, "the solver calls this function (Tfluid function) after every successful integration step". But I did not understand if it is done automatically or if I need to call it some away. Moreover:

1. how do I get the output of this Tfluid function?
2. if it is called automatically, how the parameters (t,y,flag) are defined? Is "t" just the time of the time step or is it the whole tspan defined for the PDEPE? What is the "y" variable? Is it pre-defined?
3. If I want to compute T, then I would need to do: T_n+1=T_n + (dT/dt)_n x dt. But how can I get dt and T_n?

I'm sorry, but I'm a bit lost. I would apreciate if you can help me.

[Torsten]

I also have no experience with OutputFcn. The only thing I can answer is

If I want to compute T, then I would need to do: T_n+1=T_n + (dT/dt)_n x dt. But how can I get dt and T_n

dt is given by your prescribed array "tspan" for the times when the solution should be provided. Defining

tspan = tstart:dt:tend 

in the program from which you call "pdepe", you will have to use this value for dt in OutputFcn, too. Also in OutputFcn, you could store T(n) as a persistent variable and use it in the next call to compute T(n+1) from your above formula.

[Rodrigo Resende]

@Torsten, I didn't manage to use the OutputFcn, but solved the problem in an other way. It is not an elegant solution, but it works. I did like that:

1. initiate and save the fluid temperature (Tf) and the integration time (ti)
2. the boundary function calls a function which loads Tf and ti. As it also receives the time "t", we can compute dt=t-ti and updateTf. ti and Tf are then saved again to be used in the next time step

Part of the code is the following:

In main .m file:

Tf=56;save Tf Tf;
ti=0;save ti ti;

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
    Tf=Tfluid( xl,ul,t );
    [ k,rho,cp,alpha,h ] = propriedades('interface0');
    pl=h*xl*(ul-Tf)+k*ul;
    ql=-k*xl;
    [ k,rho,cp,alpha,h ] = propriedades('interface1');
    pr=h*xr*(ur-3)+k*ur;
    qr=k*xr;
end

function [ Tf ] = Tfluid( xl,ul,t )
    [ k,rho,cp,alpha,h ] = propriedades( 'interface0' );
    load Tf
    load ti
    dt=t-ti;
    if dt>0
        dTf=-2*h*(Tf-ul)*dt/(rho*cp*xl);
        Tf=Tf+dTf
        save Tf Tf
        ti=t;
        save ti ti
    end
end

[Torsten]

If it works, ok. But keep in mind that "heatbc" - in contrast to "OutputFcn" - is not (only) called at times where the solutions ul and ur are already established, but also for "test" times and values for ul and ur in between. Despite your if-statement (if dt> 0), I have my doubts about the reliability of the results.

Further keep in mind that usually in pdepe's problem function,

f = k*DuDr 

is set.

If this is the case also in your code, don't multiply ql and qr again with k as you do it in "heatbc". Otherwise, your boundary condition will look something like

k^2*du/dr = ...

[Rodrigo Resende]

I think you're right. Maybe the results are not trustable. Unfortunately I don't know how to use OutputFcn. I had to limit MaxStep to achieve a not absurd result. Still, I'm not sure it is correct.

[Torsten]

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):

https://de.mathworks.com/matlabcentral/fileexchange/97437-pde1dm?tab=reviews

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.

Answer 3

[Torsten]

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

[Sanjay]

[T,Y] = ode45(@(t,y) odefun(t,y,nr,r,rm,wallConductivity,wallDensity,wallSpecHeat,rhoAir,cpAir,hAir,hHotAir,tempHotAir),t,y0);

i am getting error in this line of code. how I can resolve

[Torsten]

You will have to include your code and the complete error message you are getting.

The problem is stiff - so there is a reason why I used ode15s instead of ode45.

Answer 4

[Rodrigo Resende]

@Bill Greene and @Torsten,

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):

1. pde1dm.m: addParameter ---> addParamValue
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)

5. decicShampine.m: addParameter ---> addParamValue

The only change I was in doubt was if I could replace "jac" by "jac(1)".

[Bill Greene]

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