Computing temperature of a fluid inside a cylinder using PDEPE
15 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Rodrigo Resende
el 10 de Ag. de 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 comentarios
Torsten
el 11 de Ag. de 2022
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
el 11 de Ag. de 2022
Yes,I was thinking of trying pde1dm on this problem (using one ODE as you suggest) if I could understand the physics. ;-)
Respuesta aceptada
Bill Greene
el 12 de Ag. de 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
Más respuestas (3)
Torsten
el 10 de Ag. de 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 comentarios
Torsten
el 12 de Ag. de 2022
Editada: Torsten
el 12 de Ag. de 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
el 12 de Ag. de 2022
Editada: Torsten
el 12 de Ag. de 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
2 comentarios
Rodrigo Resende
el 15 de Ag. de 2022
Editada: Rodrigo Resende
el 15 de Ag. de 2022
1 comentario
Bill Greene
el 15 de Ag. de 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
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!