Computing temperature of a fluid inside a cylinder using PDEPE

Question

Rodrigo Resende 2022-8-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1776920-computing-temperature-of-a-fluid-inside-a-cylinder-using-pdepe

移动：Torsten 2024-7-12

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 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2022-8-11

@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 2022-8-11

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 2022-8-12

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1776920-computing-temperature-of-a-fluid-inside-a-cylinder-using-pdepe#answer_1025635

在 MATLAB Online 中打开

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Rodrigo Resende 2022-8-12

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 2022-8-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1776920-computing-temperature-of-a-fluid-inside-a-cylinder-using-pdepe#answer_1024285

在 MATLAB Online 中打开

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 个评论
显示 7更早的评论隐藏 7更早的评论

Rodrigo Resende 2022-8-11

编辑：Rodrigo Resende 2022-8-11

在 MATLAB Online 中打开

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

initiate and save the fluid temperature (Tf) and the integration time (ti)
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

Rodrigo Resende 2022-8-11

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 2022-8-12

编辑：Torsten 2022-8-12

在 MATLAB Online 中打开

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 2022-8-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1776920-computing-temperature-of-a-fluid-inside-a-cylinder-using-pdepe#answer_1025795

编辑：Torsten 2022-8-12

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Sanjay 2024-7-12

移动：Torsten 2024-7-12

[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 2024-7-12

移动：Torsten 2024-7-12

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 2022-8-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1776920-computing-temperature-of-a-fluid-inside-a-cylinder-using-pdepe#answer_1026500

编辑：Rodrigo Resende 2022-8-15

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

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Bill Greene 2022-8-15

在 MATLAB Online 中打开

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

请先登录，再进行评论。

Computing temperature of a fluid inside a cylinder using PDEPE

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

更多回答（3 个）

9 个评论
显示 7更早的评论隐藏 7更早的评论

2 个评论
显示无隐藏无

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

Computing temperature of a fluid inside a cylinder using PDEPE

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

更多回答（3 个）

9 个评论 显示 7更早的评论隐藏 7更早的评论

2 个评论 显示 无隐藏 无

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论

2 个评论
显示无隐藏无

1 个评论
显示 -1更早的评论隐藏 -1更早的评论