How can we define the initial conditions u0 at time t+1 as a vector that is equal to the previous initial conditions u0 at time t?

1 visualización (últimos 30 días)
Hello community,
I am currently working on a co-simulation problem that uses PDEPE toolbox to solve heat conduction.
I am utilizing EnergyPlus to simulate a cold room for a 24-hour period = endtime. The model calculates the internal temperature at each timestep, which is set to 60 seconds. This internal temperature is then employed to determine the product's temperature every minute using a 1D heat transfer approach with Robin boundary conditions (convective heat flux). Subsequently, this information is used to calculate the thermal load of the product, which is then transmitted to EnergyPlus for further analysis and processing.
And since PDEPE tool needs time span of integration, specified as a vector specifying the points at which a solution is requested, I created an array tspan_pde = 0:60:endtime to match the EnergyPlus simulation runtime.
The problem is at each EnergyPlus timestep the PDEPE tool is called and solves the conduction problem for the whole tspan_pde.
My idea is instead of solving the PDE equation for the whole tspan_pde = 0:60:endtime, I should solve it only for [0s 30s 60s] :
Tprod = Tprod0;
(1) u0 = Tprod; Q = h*s*(Tprod-Tw); will be send to EnergyPlus;
(2) EnergyPlus calculates the internal temperature Tw;
(3) Tw is used to calculate Tprod by solving PDEPE for tspan_pde = [0s 30s 60s] ==> u = [u(0) u(30) u(60)]; store Tprod = u(60)
(4) Update the initial condition u0 = u(60); Q = h*s*(Tprod-Tw); will be send to EnergyPlus;
(5) Increment t; go to (2);
So that is why I need to assign the initial conditions u0 at time t+1 to be equal to the temperature u at time t.
global rho Cp k
rho = 1000;
Cp = 888;
k = 20;
L = 1;
x = linspace(0,L,10);
endt = 4*3600;
t = linspace(0,endt,20);
m = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);
figure
plot(t,sol(:,2),t,sol(:,5),t,sol(:,8))
function [c,f,s] = heatpde(x,t,u,dudx)
global rho Cp k
c = rho*Cp;
f = k*dudx;
s = 0;
end
function u0 = heatic(x)
u0 = 20;
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur - 1;
qr = 0;
end

Respuesta aceptada

Torsten
Torsten el 29 de Mayo de 2023
Movida: Torsten el 29 de Mayo de 2023
In order to proceed, I need to assign the initial conditions u0 at time t+1 to be equal to the previous initial conditions u0 at time t.
You won't proceed this way because the solution will always remain the same as at time t=0.
  4 comentarios
Torsten
Torsten el 30 de Mayo de 2023
Something like this (where sol_old is the pdepe solution of the last 1-minute timestep) ?
function [sol_new] = pdepe_application(T_wall,sol_old)
rho = 1000;
Cp = 888;
k = 20;
L = 1;
x = linspace(0,L,10);
t = linspace(0,60,3);
heatic = @(xp)interp1(x(:),sol_old(end,:,1),xp);
m = 0;
sol_new = pdepe(m,@heatpde,heatic,@(xl,ul,xr,ur,t)heatbc(xl,ul,xr,ur,t,T_wall),x,t);
function [c,f,s] = heatpde(x,t,u,dudx)
c = rho*Cp;
f = k*dudx;
s = 0;
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
pl = ...;
ql = ...;
pr = ...;
qr = ...;
end
end
Zakaria OUAOUJA
Zakaria OUAOUJA el 8 de Jun. de 2023
Thank you very much for you answer it did really help.
clear all; clc; close all
global rho Cp k h1 h2 Tair xvec Tprodi L
rho = 1000; Cp = 8.88; k = 2; L = 1; xnodes = 5;
xvec = linspace(0,L,xnodes);
endt = 0.5*3600;
Dtmat = 60;
t = linspace(0, Dtmat, 3);
%t = linspace(0,60,1+endt/60);
m = 0; % the symmetry of the problem (slab = 0, cylindrical = 1, or spherical = 2)
h1 = 30; h2 = 35;
Tair = -20;
%Tprodi = [10 5 3 0 -5 ];
Tprodi = 25*ones(1,xnodes);
a=0; % while loop parameter
iLog=1;
while a<endt
%pde_IC = @(x)interp1(xvec,T0,x);
sol = pdepe(m,@pde_eq,@pde_IC,@pde_BC,xvec,t);
%Tprod(iLog+1,:) = [sol(3,:)]
Tprod(iLog,:) = Tprodi;
Tprodi = sol(3,:);
datat(iLog) = [a];
a = a + 60;
iLog = iLog + 1;
end
figure
plot(datat/60,Tprod(1:iLog-1,:))
xlabel('[min]')
ylabel('[°C]')
legend('node 1','node 2','node 3','node 4','node 5')
function [c,f,s] = pde_eq(x,t,u,dudx)
global rho Cp k
c = rho*Cp;
f = k*dudx;
s = 0;
end
function u0 = pde_IC(x)
global xvec Tprodi
u0 = interp1(xvec,Tprodi,x);
end
function [pl,ql,pr,qr] = pde_BC(xl,ul,xr,ur,t)
global h1 h2 Tair
pl = -h1*(ul-Tair);
ql = 1;
pr = h2*(ur-Tair);
qr = 1;
end

Iniciar sesión para comentar.

Más respuestas (0)

Productos


Versión

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by