Медленная интеграция системы ODE при определенных условиях c - PullRequest
Новая
       106

Медленная интеграция системы ODE при определенных условиях c

0 голосов
/ 06 марта 2020

Я пытаюсь решить систему ODE ниже. Я применяю конечные различия к z, поэтому я получаю систему ODE, которую я могу решить с помощью решателя, такого как ode45.

ODE system

Проблема в том, что в условиях, которые меня интересуют, интеграция идет слишком медленно (шаг по времени очень мал), и я не получаю результата. Я хотел бы знать, что является основной проблемой и есть ли способ ее решить. Я пытался:

  • с использованием различных решателей (пробовал все из них)
  • с использованием безразмерных переменных
  • с использованием функции pdepe Matlab для непосредственного решения системы PDE без применяя конечные различия сам. Возникает та же проблема.

Полный код (с использованием параметров, которые не работают) 

function transport_model()
global parameter L Di epsb Q Cfeed K tfinal npz npt F ui h

parameter = [0.044511*432.75 432.75];    % isotherm parameters
L = 0.4;            % m, column length
epsb = 0.367;       % column bulk porosity
ui = 4e4;        % m/s, velocity
Cfeed = 200;     % feed concentration
K = 0.0782          % s-1, mass transfer coefficient
tfinal = 400;      % min, final time for calculation
npz = 50;           % number of discretization points in z
npt = 20;

F = (1-epsb)/epsb;
tspan = 0:tfinal/(npt-1):tfinal;
y0 = zeros(2*npz,1);
h = L/(npz-1);

sol = ode15s(@sedo, tspan, y0);
t = sol.x;
C1 = sol.y(1:npz,:);
q1 = sol.y((npz+1):end,:);

plot(t, C1(end,:))
xlim([0 tfinal])
ylim([0 Cfeed])



function DyDt = sedo(t,y)
global Cfeed K  npz F ui h

N = npz;
y(1) = Cfeed;
DyDt = zeros(2*N,1);

% Forward finite differences
DyDt(1) = -ui * 1/(2*h)*(-y(3)+4*y(2)-3*y(1)) - F * K * (y(1) - ilangmuir(y(N+1)));
DyDt(N+1) = K * (y(1) - ilangmuir(y(N+1)));

% Central finite differences
for i=2:N-1
    DyDt(i) = -ui * 1/(2*h)*(y(i+1)-y(i-1)) - F * K * (y(i) - ilangmuir(y(N+i)));
    DyDt(N+i) = K * (y(i) - ilangmuir(y(N+i)));
end

% Backward finite differences
DyDt(N) = -ui * 1/(2*h)*(3*y(N)-4*y(N-1)+y(N-2)) - F * K * (y(N) - ilangmuir(y(N+N)));
DyDt(2*N) = K * (y(N) - ilangmuir(y(N+N)));



function c = ilangmuir(q)
% langmuir solved for c
global parameter

a = parameter(1);
b = parameter(2);

c = q /(a - b * q);

Однако, если я пытаюсь смоделировать другой случай ( разные условия) я получаю правильный ответ. Поэтому я предполагаю, что это численная проблема, вызванная условиями, которые я пытаюсь смоделировать, а не ошибка в модели. Например, если я использую условия ниже, я получаю правильный ответ. 

parameter = [0.044511*432.75 432.75];    % isotherm parameters
L = 0.4;            % m, column length
epsb = 0.367;       % column bulk porosity
ui = 0.0056;        % m/s, velocity
Cfeed = 0.0025;     % feed concentration
K = 0.0782          % s-1, mass transfer coefficient
tfinal = 4000;      % min, final time for calculation
npz = 50;           % number of discretization points in z
npt = 20;

enter image description here

...