Solving PDE and ODE coupled system with varying boundary conditions?

I'm trying to recreate a model of desiccant wheel from this paper using matlab:

https://www.researchgate.net/publication/241826275_Performance_Modeling_of_Desiccant_Wheels_1_Model_Development

It simplifies the desiccant wheel into one single 1D channel that goes through periods of process and regeneration.

The balance equations can then be written as:

which should be an system of pdes (eq.1 and eq.3) and odes (eq.2 and eq.4 if we set k_sub = 0).

The only variables that changes with x and time here are , with calculated from and using the isotherm function.

The boundary varies overtime with being set at either x=0 or x=L depending if the channel is in process or regeneration region.

I've did some research and found that pdede cannot solve a system containing odes so I attempted to use the pde1dm library to code it. Here is my current code:

function DW_1D
    nx = 50;
    L = 0.20;
    x = linspace(0,L,nx);
    t = linspace(0,100,50);
    xOde = x;

    m = 0;
    
    [sol,odeSol] = pde1dm(m,@pdefun,@pdeic,@pdebc,x,t,@odefun,@odeic,xOde);

    rho_vg = sol(:,:,1);
    t_g = sol(:,:,2);
    
    s = surf(x,t,rho_vg);
    xlabel('x');
    ylabel('t');
    s.EdgeColor = 'none';
end


function f=odefun(t,v,vdot,x,u,DuDx)
    rho_vg = u(1);
    t_g = u(2);
    gamma_m = v(1);
    t_m = v(2);
    rho_vm = isotherm(gamma_m, t_m);

    h = 43.3;
    h_m = 0.036;

    flow_speed = 1.3;
    channel_a = 0.0015;
    channel_b = 0.0034;
    P = (channel_a + channel_b)*2;
    A = channel_a * channel_b;
    m_thickness = 65E-6;
    sub_thickness = 75E-6;
    A_m = P * m_thickness;
    A_sub = P * sub_thickness;
    Cp_m = 1000;
    Cp_sub = 900;
    delta_abs = 2791000;
    rho_m = 700;
    rho_sub = 500;

    rho_g = 1.293;
    Cp_g = 1006;

    f = [h_m*P*(rho_vg-rho_vm)-A_m*rho_m*vdot(1); h_m*P*(rho_vg-rho_vm)*delta_abs - h*P*(t_m-t_g) - (rho_m*A_m*Cp_m + rho_sub*A_sub*Cp_sub)*vdot(2)];
end

function v0 = odeic(t0)
    v0 = [0.1; 300];
end

function [c,f,s] = pdefun(x,t,u,dudx,v,vdot) % Equation to solve
    rho_vg = u(1);
    t_g = u(2);
    gamma_m = v(1);
    t_m = v(2);
    rho_vm = isotherm(gamma_m, t_m);

    h = 43.3;
    h_m = 0.036;

    flow_speed = 1.3;
    channel_a = 0.0015;
    channel_b = 0.0034;
    P = (channel_a + channel_b)*2;
    A = channel_a * channel_b;
    m_thickness = 65E-6;
    sub_thickness = 75E-6;
    A_m = P * m_thickness;
    A_sub = P * sub_thickness;
    Cp_m = 1000;
    Cp_sub = 900;
    delta_abs = 2791000;
    rho_m = 700;
    rho_sub = 500;

    rho_g = 1.293;
    Cp_g = 1006;

    c = [1; 1];
    f = [0; 0];
    s = [-flow_speed*dudx(1)+(-(h_m*P)/A)*(rho_vg-rho_vm);
        -flow_speed*dudx(2)+((h*P)/(A*rho_g*Cp_g))*(t_m-t_g)];   
end
% ---------------------------------------------
function u0 = pdeic(x) % Initial Conditions
    u0 = [0.005; 300];
end
% ---------------------------------------------
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions
    pl = [ul(1)-0.005; ul(2)-340];
    ql = [0; 0];
    pr = [0; 0];
    qr = [1; 1];
end
% ---------------------------------------------
function rho_vm = isotherm(gamma_m, Tm) % Isotherm
    gamma_max = 0.36;
    c = 1;
    relative_humidity = (gamma_m*c)/(gamma_max-gamma_m+gamma_m*c);
    rho_vm = relative_humidity/(4.09E-9 * Tm * exp(5196/Tm));
end

My first problem is that the function seems to run into an error when I use only two boundary conditions (setting to 340 and 0.005 respectively) such as the following, without setting qr to [1;1] matlab will giveout an error "Failure of consistent initialization algorithm: There are 4 algebraic equations but the rank is only 2.". When I set qr to [1;1] it seems to run fine.

function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions
    pl = [ul(1)-0.005; ul(2)-340];
    ql = [0; 0];
    pr = [0; 0];
    qr = [0; 0];
end

My second problem is that I can't seem to figure out how to change the boundary condition while retaining the current state. I tried to extract the variables so I can rerun the solver using the last state with different boundary condition but the solver only outputs the solved Pde solution as an array at different x but not the Ode solution. The ode solution only contains one value when the value should vary at different x position.

Can anyone help me with these two problems?

Thanks!

Edit: The paper also specifies a set of boundary conditions for (calculated from and ) and but I have no idea how to write it:

[Torsten]

Edit: The paper also specifies a set of boundary conditions for (calculated from and ) and but I have no idea how to write it.

A boundary condition of the above type for t_m is possible if you include the second-order spatial derivative d^2t_m/dx^2.

Setting a boundary condition for rho_vm is not possible because it is a quantity deduced from gamma_m and t_m.

[Andy Fang]

Thank you so much!

Answer 1

[Torsten]

If you want to change a boundary condition, use a second function "fun" and call ode15s with the solution obtained so far and this new function.

nx = 50;
L = 0.20;
x = linspace(0,L,nx).';
t = linspace(0,100,50);
u0 = [0.005; 300]; %rho_vg,t_g
v0 = [0.1; 300];  %gamma_m,t_m
y0 = [u0(1)*ones(nx,1);340;u0(2)*ones(nx-1,1);v0(1)*ones(nx,1);v0(2)*ones(nx,1)];
[T,Y] = ode15s(@(t,y)fun(t,y,x,nx),t,y0);
figure(1)
plot(x,Y(end,1:nx))

figure(2)
plot(x,Y(end,nx+1:2*nx))

figure(3)
plot(x,Y(end,2*nx+1:3*nx))

figure(4)
plot(x,Y(end,3*nx+1:4*nx))

function dydt = fun(t,y,x,nx)
  rho_vg = y(1:nx);
  t_g = y(nx+1:2*nx);
  gamma_m = y(2*nx+1:3*nx);
  t_m = y(3*nx+1:4*nx);
  rho_vm = isotherm(gamma_m,t_m);

  h = 43.3;
  h_m = 0.036;
  flow_speed = 1.3;
  channel_a = 0.0015;
  channel_b = 0.0034;
  P = (channel_a + channel_b)*2;
  A = channel_a * channel_b;
  m_thickness = 65E-6;
  sub_thickness = 75E-6;
  A_m = P * m_thickness;
  A_sub = P * sub_thickness;
  Cp_m = 1000;
  Cp_sub = 900;
  delta_abs = 2791000;
  rho_m = 700;
  rho_sub = 500;
  rho_g = 1.293;
  Cp_g = 1006;

  k_sub = 0.4;   % Setting to close the equations
  
  d_rho_vg_dt = zeros(nx,1);
  d_t_g_dt = zeros(nx,1);
  d_gamma_m_dt = zeros(nx,1);
  d_t_m_dt = zeros(nx,1);
  
  %1st order
  %d_rho_vg_dt(1) = 0.0;
  %d_rho_vg_dt(2:nx) = -flow_speed*(rho_vg(2:nx)-rho_vg(1:nx-1))./(x(2:nx)-x(1:nx-1))+...
  %                    (-(h_m*P)/A)*(rho_vg(2:nx)-rho_vm(2:nx));
  %d_t_g_dt(1) = 0.0;
  %d_t_g_dt(2:nx) = -flow_speed*(t_g(2:nx)-t_g(1:nx-1))./(x(2:nx)-x(1:nx-1))+...
  %                 ((h*P)/(A*rho_g*Cp_g))*(t_m(2:nx)-t_g(2:nx));

  %2nd order
  d_rho_vg_dt(1) = 0.0;
  d_rho_vg_dt(2:nx-1) = -flow_speed*(rho_vg(3:nx)-rho_vg(1:nx-2))./(x(3:nx)-x(1:nx-2))+...
                        (-(h_m*P)/A)*(rho_vg(2:nx-1)-rho_vm(2:nx-1));
  d_rho_vg_dt(nx) = -flow_speed*(rho_vg(nx)-rho_vg(nx-1))./(x(nx)-x(nx-1))+...
                    (-(h_m*P)/A)*(rho_vg(nx)-rho_vm(nx));
  
  d_t_g_dt(1) = 0.0;
  d_t_g_dt(2:nx-1) = -flow_speed*(t_g(3:nx)-t_g(1:nx-2))./(x(3:nx)-x(1:nx-2))+...
                     ((h*P)/(A*rho_g*Cp_g))*(t_m(2:nx-1)-t_g(2:nx-1));
  d_t_g_dt(nx) = -flow_speed*(t_g(nx)-t_g(nx-1))./(x(nx)-x(nx-1))+...
                 ((h*P)/(A*rho_g*Cp_g))*(t_m(nx)-t_g(nx));
  
  d_gamma_m_dt = h_m*P*(rho_vg-rho_vm)/(A_m*rho_m);
  
  % Include d^2t_m/dx^2 term with d_tm/dx = 0 at x = 0 and x = L
  d_t_m_dt(1) = (k_sub*A_sub*(t_m(2)-t_m(1))/(x(2)-x(1))/((x(1)+x(2))/2-x(1))+...
               h_m*P*(rho_vg(1)-rho_vm(1))*delta_abs - h*P*(t_m(1)-t_g(1)))/...
               (rho_m*A_m*Cp_m + rho_sub*A_sub*Cp_sub);
  d_t_m_dt(2:nx-1) = (k_sub*A_sub*((t_m(3:nx)-t_m(2:nx-1))./(x(3:nx)-x(2:nx-1))-...
                      (t_m(2:nx-1)-t_m(1:nx-2))./(x(2:nx-1)-x(1:nx-2)))./...
                      ((x(3:nx)+x(2:nx-1))/2 - (x(2:nx-1)+x(1:nx-2))/2) +...
                     h_m*P*(rho_vg(2:nx-1)-rho_vm(2:nx-1))*delta_abs - h*P*(t_m(2:nx-1)-t_g(2:nx-1)))/...
                     (rho_m*A_m*Cp_m + rho_sub*A_sub*Cp_sub);
  d_t_m_dt(nx) = (k_sub*A_sub*(-(t_m(nx)-t_m(nx-1))/(x(nx)-x(nx-1))/(x(nx)-(x(nx)+x(nx-1))/2))+...
                  h_m*P*(rho_vg(nx)-rho_vm(nx))*delta_abs - h*P*(t_m(nx)-t_g(nx)))/...
                  (rho_m*A_m*Cp_m + rho_sub*A_sub*Cp_sub);
  
  %d_t_m_dt = (h_m*P*(rho_vg-rho_vm)*delta_abs - h*P*(t_m-t_g))/(rho_m*A_m*Cp_m + rho_sub*A_sub*Cp_sub);
  dydt = [d_rho_vg_dt;d_t_g_dt;d_gamma_m_dt;d_t_m_dt];
  
end

% ---------------------------------------------
function rho_vm = isotherm(gamma_m, Tm) % Isotherm
    gamma_max = 0.36;
    c = 1;
    relative_humidity = (gamma_m*c)./(gamma_max-gamma_m+gamma_m*c);
    rho_vm = relative_humidity./(4.09E-9 * Tm .* exp(5196./Tm));
end