Clamping solution to a BVP solver

Hey all,

below you can find a standard solution to solve steady-state Fickian transport with bvp4c. I am trying in the moment to impose a c_max treshhold to my function (or boundary condtion) so the solution for c is not exceeding the preset c_max. I do not want to impose a second boundary to the c itself.

The solution should limit c without knowing that c_sat is a threshold for c through the boundary condition.

Has anybody done something like this before or could give me a hint on how this clamp could be done using bvp4c?

main_clamp_boundary()

function main_clamp_boundary()
clear all; close all; clc;
    % Parameters
    c_start   = 1e-3;            % base concentration
    c_max     = 50*c_start;      % threshold
    D         = 1e-12;
    S         = 1e-2;               % (positive => 'production' or 'source')
    L         = 5e-6;
    k_penalty = 1e9;             % large penalty factor

    % Mesh and initial guess
    xmesh   = linspace(0, L, 100);
    solinit = bvpinit(xmesh, @guess);

    % Solve using bvp4c
    sol = bvp4c(@odefun, @bcfun, solinit);

    % Extract solution
    x    = sol.x;
    c    = sol.y(1,:);
    dcdx = sol.y(2,:);

    figure()
    plot(x, c, 'LineWidth', 2)
    hold on
    yline(c_max, '--r', 'c_{max}', 'LineWidth', 1)
    xlabel('x'), ylabel('c(x)')
    title('Concentration with Soft Clamp at x=L (Boundary Condition)')
    legend('c(x)','c_{max}','Location','best')

    figure()
    plot(x, dcdx, 'LineWidth', 2)
    xlabel('x'), ylabel('dc/dx')
    title('Concentration Gradient')
    grid on

    function dydx = odefun(x,y)
 
        c    = y(1);
        dcdx = y(2);

        dydx = zeros(2,1);
        dydx(1) = dcdx;    
        dydx(2) = S / D;     
    end

    function res = bcfun(ya, yb)

       res = [ya(1);
              ya(2)];
    end

%%  Penalty factor

    % function dydx = odefun(x,y)
    %     % ODE system: y(1) = c, y(2) = dc/dx
    %     c    = y(1);
    %     dcdx = y(2);
    %     k_penatly = 9e9;
    % 
    %     penalty = k_penalty * max(0, c - c_max);
    %     dydx(2) = S / D - k_penalty;     % c'' = S/D (constant source)
    % end

    function g = guess(x)
        g = [c_start; 0];
    end
end

[Torsten]

Do you mean a mixed boundary condition like -D*dc/dn = beta*(c-c_max) with n being the outer normal ?

[Thanh Hoang]

Hi Torsten,

I was playing with the code and did not copied in the original problem. The graphs now are representing the initial state which I want to solve.

The concentration is exceeding the c_max that I set and I just try out ideas how I could limit this without applying a second Drichlet boundary condition to c.

A Robin type boundary condition should be actually already applied to the problem, however I want to clamp c on the other side of the domain. I am not sure if an additional boundary condition could be applied here, though my question was if this clamping would also be possible without setting a boundary.

[Torsten]

Your equation has solution c(x) = S/(2*D)*x^2. So if you want to limit c without changing a boundary condition, you can either change S or D or L.