Pressure Swing Adsorption: Solve coupled partial differential equation system with ODE-Solver and Method of Lines

Dear all,

I want to solve the following PDE, which describes Fluid Flow in space and time:

I want to apply the "Method of Lines" and ODE-Solver in Matlab to solve the PDE, in order to extend the PDE to coupled PDE-Systems shown in the appendix ("Equation system with BC 2" & "Equation System with BC 3"). My corresponding MATLAB-code to solve the simplest fluid flow is:

 % clear all; clc;
%Parameters: 
v=0.1; %velocity m/s
%Boundary Conditions:  
C_BC=5;% BC for space x=0
L=1; %space length 
n=50; %space points
dx=L/n; %space intervall  
%Initial Condition: 
for i=1:n 
     u0(i)=0; 
end
t0=0; %lower bound for time 
tf=8; %upper bound for time 
m=50; %time points 
tspan=linspace(t0,tf,m); %time vector  
[t,u]=ode15s(@(t,u)Fluid_Flow(u,n,v,C_BC,dx),tspan,u0);
for i=1:n 
    if i==1 
        space(i) = 0; 
    else 
        space(i) = space(i-1)+dx ; 
    end 
end 
%concatenate the solution variable/vector 
for i=1:n 
    if i==1 
    C=u(:,1) ; 
    else 
    C=cat(2,C,u(:,i));
    end 
end 
h=surf(space,t,C);
set(h,'LineStyle','none');
xlabel('space x')
ylabel('time t')
zlabel('Concentration C')
function ut=Fluid_Flow(u,n,v,C_BC,dx)
    for i=1:n     
        if(i==1)  
            ut(i)= -1*v*(u(i)-C_BC)/dx;
        elseif (i==n) 
            ut(i)=0;
        else
            ut(i)= -1*v*(u(i)-u(i-1))/dx;
        end 
    end
    ut=ut';
end

I expect to get a "sharp" plug-flow of fluid, since there is no diffusion-term in the pde. But my visualization shows a "dispersed" flow as you can see in the red marked area in the following graphic:

I don't know how to simulate the fluid-flow correctly. Maybe the boundary conditions doesn't fit or there is another profound failure. I would appreciate to get some suggestions and information, e.g. how to set the BC correctly.

I want to solve the coupled PDE in principle by extending the vector u(i) within the ODE-Solver like in the following function (random function without meaning):

    % function ut=Fluid_Flow(u,n,v,C_BC,dx)
    for i=1:n     
        if(i==1)  
            ut(i)   = -1*v*(u(i)-C_BC)/dx;
            ut(n+i) = -v*(u(n+i)+C_BC)/dx;
        elseif (i==n) 
            ut(i)   =0;
            ut(n+i) =0;
        else
            ut(i)   = -1*v*(u(i)-u(i-1))/dx;
            ut(n+i) = -v*(u(i)-u(i-1))/(dx^2); 
        end 
    end
    ut=ut';
    end

Initial condition have to be extended in the main-programm as well:

   %Initial Condition: 
for i=1:n 
     u0(i)   = 0; 
     u0(n+i) = 0;
end

I'm not sure, if the coupled PDE-system works, if the coupling is executed this way. I would appreciate, if someone could confirm it.

Thank You.

Answer 1

[Torsten]

You need a second-order discretization in space to avoid these artificial diffusion effects.

Take a look at

http://www.nag.co.uk/numeric/mb/nagdoc_mb/manual_25_1/pdf/d03/d03intro.pdf

chapter 3.7 for adequate solvers.

Best wishes

Torsten.

[Robert Anhalt]

Thanks a lot for Your suggestion!!

I've reviewed Your provided link! I could find further helpful literature, which deals with artificial diffusion respectively numerical diffusion.

Artificial diffusion is limited, when only the number of space intervalls is increased as you can see in the graphic (I set n=15000 in my code):

But the computation for this method is time-consuming. Hence I will implement a more efficient upwind scheme for space discretization. I hope this will work.

Kind regards

Robert

[Torsten]

By the way:

You don't need to distinguish between 1<i<n and i=n in your discretization:

function ut=Fluid_Flow(u,n,v,C_BC,dx)
    for i=1:n     
        if(i==1)  
            ut(i)= -1*v*(u(i)-C_BC)/dx;
        else
            ut(i)= -1*v*(u(i)-u(i-1))/dx;
        end 
    end
    ut=ut';
end

Best wishes

Torsten.