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?

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

Answer 1

[Torsten]

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.

[Zakaria OUAOUJA]

Thank you for your comment

I've updated the description to make it clearer

[Torsten]

It's still unclear (at least for me).

What information do you get from EnergyPlus that is needed for pdepe or what information do you get from pdepe that is needed in EnergyPlus ?

And you want to run EnergyPlus for 1 minute, give information to pdepe from EnergyPlus and run it for one minute, again give information to EnergyPlus from pdepe and run it for one minute and so on ?

[Torsten]

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]

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