Time interation error while trying to solve a partial differential equation which describes pore diffusion resistance combined with with surface kinetics.

I have attached the additional information (problem statement 2). The governing equation and the boundary conditions are given in it. I used pdepe function to solve the PDE. Here is the code :

function concentration_profiles_time()
    % Define parameters
    Rp = 1e-3; % Radius of the pellet (m)
    epsilon_p = 0.5; % Porosity
    De = 1.5e-10; % Effective diffusivity (m^2/s)
    kf = 2.8e-3; % Mass transfer coefficient (m/s)
    kR = 5.0e-2; % Reaction rate constant (1/s)
    Cb = 1; % Bulk concentration (arbitrary units)
    
    % Define spatial and temporal discretization
    r = linspace(0, Rp, 50); % Radial discretization
    t = linspace(0, 100, 100); % Time discretization (s)

    % Solve PDE using pdepe with options for integration tolerance
    m = 2; % Spherical symmetry
    sol = pdepe(m, @pde_eqn, @initial_condition, @boundary_conditions, r, t);
    
    % Plot results (check the number of steps in the solution)
    figure;
    num_steps = size(sol, 1);  % Number of time steps in the solution
    for i = 1:min(10, num_steps) % Ensure we don't exceed the available number of steps
        plot(r, sol(i, :), 'DisplayName', sprintf('t = %.1f s', t(i)));
        hold on;
    end
    xlabel('Radius (m)');
    ylabel('Concentration (C_p)');
    legend;
    title('Concentration Profiles over Time');
    grid on;
    hold off;
    
    % Nested functions
    function [c, f, s] = pde_eqn(r, t, Cp, dCdr)
        % PDE coefficients
        c = epsilon_p; % Time-dependent term coefficient
        f = De * dCdr; % Diffusion term coefficient
        s = -kR * Cp; % Source term
    end

    function Cp0 = initial_condition(r)
        % Initial condition
        Cp0 = zeros(size(r)); % Initial concentration is zero everywhere
    end

    function [pl, ql, pr, qr] = boundary_conditions(xl, ul, xr, ur, t)
        % Boundary conditions
        % At r = 0 (center of the sphere), Neumann BC: dCp/dr = 0
        pl = 0;
        ql = 1;
        
        % At r = Rp (surface of the sphere)
        pr = -kf * (Cb - ur) - kR * ur;
        qr = De;
    end
end

This is the message I receive when running the code.

Warning: Failure at t=1.090867e-12. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.231174e-27) at time t.

> In ode15s (line 625)

In pdepe (line 290)

In Testing (line 16)

Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00.

> In pdepe (line 306)

In Testing (line 16)

It would would be really helpful if someone could tell me what I'm doing wrong here

Answer 1

[Torsten]

One error, one problem:

The error:

qr = De;

must be replaced by

qr = 1;

because the boundary condition is given by p+q*f = 0 where f = De*dCdr in your case. Thus De is already included in your f.

The problem:

Your reaction rate constant is chosen too large. The substance produced at the boundary at r = R_p cannot be transported to the interior of the pellet fast enough because of the small diffusion coefficient. That's why the concentration at the boundary grows to infinity (see below). That's at least what I'm guessing.

concentration_profiles_time()

function concentration_profiles_time()
    % Define parameters
    Rp = 1e-3; % Radius of the pellet (m)
    epsilon_p = 0.5; % Porosity
    De = 1.5e-10; % Effective diffusivity (m^2/s)
    kf = 2.8e-3; % Mass transfer coefficient (m/s)
    kR = 5.0e-2; % Reaction rate constant (1/s)
    Cb = 1; % Bulk concentration (arbitrary units)
    
    % Define spatial and temporal discretization
    r = linspace(0, Rp, 50); % Radial discretization
    t = linspace(0, 10, 11); % Time discretization (s)

    % Solve PDE using pdepe with options for integration tolerance
    m = 2; % Spherical symmetry
    sol = pdepe(m, @pde_eqn, @initial_condition, @boundary_conditions, r, t,odeset('RelTol',1e-10,'AbsTol',1e-10));
    
    % Plot results (check the number of steps in the solution)
    figure;
    num_steps = size(sol, 1);  % Number of time steps in the solution
    for i = 1:num_steps % Ensure we don't exceed the available number of steps
        plot(r, sol(i, :), 'DisplayName', sprintf('t = %.1f s', t(i)));
        hold on;
    end
    xlabel('Radius (m)');
    ylabel('Concentration (C_p)');
    legend;
    title('Concentration Profiles over Time');
    grid on;
    hold off;
    
    % Nested functions
    function [c, f, s] = pde_eqn(r, t, Cp, dCdr)
        % PDE coefficients
        c = epsilon_p; % Time-dependent term coefficient
        f = De * dCdr; % Diffusion term coefficient
        s = -kR * Cp; % Source term
    end

    function Cp0 = initial_condition(r)
        % Initial condition
        Cp0 = zeros(size(r)); % Initial concentration is zero everywhere
    end

    function [pl, ql, pr, qr] = boundary_conditions(xl, ul, xr, ur, t)
        % Boundary conditions
        % At r = 0 (center of the sphere), Neumann BC: dCp/dr = 0
        pl = 0;
        ql = 1;
        
        % At r = Rp (surface of the sphere)
        pr = -kf * (Cb - ur) - kR * ur*5.5e-2;
        qr = 1;
    end
end

[Aditya]

Thanks a lot for your guidance