Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab Pdepe solver for diffusion equation in a packed column

  1. Mar 30, 2016 #1
    Consider I have a packed column of length L filled with known characteristic adsorbent. I am putting a mixture of N components in it and I am solving for concentration of each component in mobile phase at the outlet of the column. The equations which are to be generalised are as follows: An example of the 2 component systems is stated in the paper below which can help you to understand the equations better if following is unclear.http://www.sciencedirect.com/science/article/pii/S1369703X10000069?np=y

    for component i:N:
    Dcoe*∂2Ci/∂x2−v*∂Ci/∂x−∂qi/∂t=∂Ci/∂t
    ∂q1∂t=ki[Ciqmi(1−Σ(qi/qmi))−kd1q1]

    [PLAIN]http://www.sciencedirect.com/sd/blank.gifwhere [Broken] the initial conditions are
    t=0,x=0,Ci(0,0)=Cfi,qi(0,0)=0,
    http://www.sciencedirect.com/sd/blank.gif
    t=0,x≠0,Ci(x,0)=0,qi(x,0)=0.

    [PLAIN]http://www.sciencedirect.com/sd/blank.gifand [Broken] the boundary conditions are
    x=0,Ci(0,t)=Cfi,qi(0,t)=0
    http://www.sciencedirect.com/sd/blank.gif
    x=L, ∂Ci/∂x=0,∂qi/∂x=0.

    I am solving these equations by using pdepe solver.

    Following 4 files are used to solve pdepe: mcl_brkthrupdemain.m (pdepe main solution function), mcl_brkthrupde (stating equations), mcl_brkthruic.m (initial conditions), mcl_brkthrubc(boundary conditions)
    %% file mcl_brkthrupdemain is used to define the input parameters and the pdepe main functions for solving the equations.
    Code (Text):

    function mcl_brkthrupdemain
    mcl_globalinc;
    m=0;
    %% Importing system parameters from an excel sheets.
    system_matrix=xlsread('mcl system physical constants.xlsx');

    % input('enter the value for velocity V in cm/min:');
    vel = system_matrix(1,1);

    % input('enter the value for Length of the column in cm:');
    L= system_matrix(1,2);

    % input('enter the value for column diameter in cm:');
    % Not used currently
    dc= system_matrix(1,3);

    % input('enter the value of porosity:');
    epsilon= system_matrix(1,4);

    %% importing parameters for component specific constants from excel sheets
    matrix_in = xlsread('mcl Excelimport.xlsx'); % Main input matrix file from user
    [nComp,nCin] = size(matrix_in); % defining number of components

    k1 = zeros(nComp,1);
    kd = zeros(nComp,1);
    qm = zeros(nComp,1);
    Cin = zeros(nComp,1);
    Dcoe = zeros(nComp,1); % input('enter the value for diffusion coefficient Dcoe in cm^2/min:');

    for w=1:nComp
    k1(w) = matrix_in(w,1);
    kd(w) = matrix_in(w,2);
    qm(w) = matrix_in(w,3);
    Cin(w) = matrix_in(w,4);
    Dcoe(w) = matrix_in(w,5);
    end
    %% Defining mesh sizes for both space and time.
    N=50;
    Run_time=800;
    % N =input('enter the value for number of meshing intervals:\n');
    % Run_time= input('enter the value for maximum time of run:\n');

    x=linspace(0,L,N); % spatial meshing
    t=linspace(0,Run_time,N); % time meshing
    tend=300;
    %% solution and plotting
    sol=pdepe(m,@mcl_brkthrupde,@mcl_brkthruic,@mcl_brkthrubc,x,t);
    colorvec= hsv(5);
    for g=1:nComp
    C=sol(:,:,2*g-1);
    solid_conc=sol(:,:,2*g);
    Cfr=C/matrix_in(g,4);
    disp(C);
    disp(solid_conc);
    % surface plots for conentrations
    %surf(x,t,C);
    %title(sprintf('surface plot for concentration of component %d',g))

    % A solution profile for conc. vs length for all components.
    figure, plot(x,C(end,:), 'color',colorvec(g,:))
    title(sprintf('solution at runtime end for Conc %d vs distance travelled',g))
    xlabel('distance (cm)')
    ylabel(sprintf ('component %d concentration',g))
    hold on

    %Plot conc. vs. time for component 1
    figure, plot(t,Cfr)
    title(sprintf('breakthrough for component %d at all x',g))
    xlabel('Time (min)')
    ylabel(sprintf('Component %d Concentration (mg/ml)',g))

    % plot end concentration with time
    figure, plot(t,Cfr(:,end))
    title(sprintf('breakthrough for Component %d',g))
    xlabel('Time (min)')
    ylabel(sprintf('Concentration of component%d (mg/ml)',g))
    end
    end

    %% file mcl_brkthrupde is used to define equations.
    Code (Text):
    function [c,f,s]= mcl_brkthrupde(x,t,u,DuDx) % inputs for equation coefficients
    %% Redefining global variables
    %% Defining system constants
    mcl_globalinc;
    neq=2*nComp;
    ddf=x;
    %% Usual variables
    conc =zeros(nComp,1);
    q =zeros(nComp,1);
    dcdx =zeros(nComp,1);
    dqdx =zeros(nComp,1);
    dqdt =zeros(nComp,1);

    for w=1:nComp
    conc(w,1) = u(2*w-1);
    q(w,1) = u(2*w);
    dcdx(w,1) = DuDx(2*w-1);
    dqdx(w,1) = DuDx(2*w);
    end
    %% writing c term and putting it into a column vector
    c=ones(neq,1);
    %disp(c);
    %% writing f term and putting it into a column vector
    f=zeros(neq,1);
    count1=1;
    while count1<(neq) % f term for all the equations
    for p=1:nComp
    f(count1,1)=Dcoe(p).*dcdx(p,1);
    f(count1+1,1)=0.*dqdx(p,1);
    count1=count1+2;
    end
    end
    disp(f);
    %% finding intermediate summation for q/qmax
    s=zeros(neq,1);
    SUM=0;
    for w=1:nComp
    SUM=SUM+q(w,1)/qm(w);
    end
    disp(SUM);
    for r=1:nComp
    dqdt(r,1)=k1(r)*conc(r,1)*qm(r)*(1-SUM)-k1(r)*kd(r).*q(r,1);
    end
    count2=1;
    while count2<(neq) % f term for all the equations
    for b=1:nComp
    s(count2,1)=-vel*dcdx(b,1)-((1-epsilon)/epsilon).*dqdt(b,1);
    s(count2+1,1)=dqdt(b,1);
    count2=count2+2;
    end
    end
    disp(s);
    end
    %% mcl_brkthruic describes initial conditions
    Code (Text):
    function u0= mcl_brkthruic(x) %function for inlet conditions
    mcl_globalinc;
    %% defining initial conditions for concentration for x=0 and 0<x<L
    u0=zeros(2*nComp,1);
    if x==0
    for i=1:nComp
    u0(2*i-1,1)= Cin(i);
    u0(2*i,1)= 0;
    end
    end
    disp(u0);
    end

    %% mcl_brkthrubc describes boundary conditions

    Code (Text):
    function [pl,ql,pr,qr]= mcl_brkthrubc(xl,ul,xr,ur,t) % function for boundary conditions
    %% calling required variable set
    mcl_globalinc;

    neq=2*nComp;
    %% defining usual variables
    pl=zeros(neq,1);
    ql=zeros(neq,1);
    pr=zeros(neq,1);
    qr=zeros(neq,1);

    cleft=zeros(nComp);
    qleft=zeros(nComp);
    cright=zeros(nComp);
    qright=zeros(nComp);

    %% finding pl (left boundary) for terms without differential term and putting into a column vector
    for i=1:nComp
    cleft(i)=ul(2*i-1);
    qleft(i)=ul(2*i);
    cright(i)= ur(2*i-1);
    qright(i)=ur(2*i);
    end
    for j=1:nComp;
    pl(2*j-1,1)=cleft(j)-Cin(j); %left for conc in liquid phase
    pl(2*j,1)=qleft(j); % left for conc in solid phase
    ql(2*j-1,1)=0; % left for differential term for conc in liquid phase
    ql(2*j,1)=0; % left for differential term for conc in solid phase

    pr(2*j-1,1)=cright(j); %right for conc in liquid phase
    pr(2*j,1)=qright(j); % right for conc in solid phase
    qr(2*j-1,1)=1; % right for differential term for conc in liquid phase
    qr(2*j,1)=1; % right for differential term for conc in solid phase
    end
    disp(pl);
    disp(ql);
    disp(pr);
    disp(qr);
    end

    As you can see that the paper mentions a bi component system and I am trying to make it a multicomponent system.
    This program is taking values from input xlsx files "mcl_system physical constants.xlsx" and "mcl Excelimport.xlsx". You can take the values you like to start with.
    It is giving me the profile but not the correct profile.
    Can anyone please have a look?
    Let me know if anything is unclear.


    Thanks for your time.
    -SimOpera
     
    Last edited by a moderator: May 7, 2017
  2. jcsd
  3. Apr 5, 2016 #2
    Thanks for the post! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post?
     
  4. Apr 24, 2016 #3
    No, Can you please pass on the query further?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Pdepe solver for diffusion equation in a packed column
  1. MatLab pdepe (Replies: 4)

Loading...