Runge Kutta 4th Order Method for an Equation

The following is an supercritical CO2 equation of an oil from leaves:

Where

• W = CO2 Mass flowrate (constant)
• ρ = solvent density (constant)
• ϵ = bed porosity (constant)
• V = extrtactor volume (constant)
• ti = internal diffusion diameter (constant)
• n = number of stages (desired to reach a stage of 10 hence n starts at 1 and ends at 10)
• Cn = fluid phase concentration in the nth stage
• Cn_bar = solid phase concentration in the nth phase
• Cn_bar* = Cn/0.2

The inital conditions are t = 0, Cn = 0 and Cn_bar = 0

I would like to know what is the best way to set up the equation and code in order for it execute the 4th order runge kutta method on Matlab so as to achieve Cn at the 10th stage. Also note the initial oil concentration in the solid phase was 9.0 kg/m3. This is the code that was developed so far, however it is not running.

clear; clc;

% Parameters
h = 20; % Step size
tfinal = 500; % Solve from t=(0,tfinal)
e = 0.38; % Bed porosity
rho = 126; % Solvent (CO2) density
W = 0.95; % CO2 flow rate, kg/h
V = 0.0004; % Extractor volume, m^3
Di = 6 * 10^-13; % Diffusion coefficient
n = 10; % Number of stages/bed divisions
ti = 1736.111; % Internal diffusion time

% Initial Condition
t(1) = 0;
c(1) = 0;
cbar(1) = 9;
c_in(1) = 0;

% Define the ODE functions
f1 = @(t,c,cbar) -1/ti * (cbar - c/0.2);
f2 = @(t,c,cbar) -(((1 - e) * V/n * f1(t,c,cbar)) + W/rho * (c - c_in))/e * V/n;

% RK4 Loop for 9 iterations
for iter = 1:10
    % RK4 Loop for solving ODEs
    for i = 1:ceil(tfinal/h)
        t(i+1) = t(i) + h;

        % Update cbar
        k1_1 = f1(t(i), c(i), cbar(i));
        k2_1 = f1(t(i) + 0.5*h, c(i) + 0.5*h*k1_1, cbar(i) + 0.5*h*k1_1);
        k3_1 = f1(t(i) + 0.5*h, c(i) + 0.5*h*k2_1, cbar(i) + 0.5*h*k2_1);
        k4_1 = f1(t(i) + h, c(i) + h*k3_1, cbar(i) + h*k3_1);
        cbar(i+1) = cbar(i) + h/6*(k1_1 + 2*k2_1 + 2*k3_1 + k4_1);

        % Update c
        k1_2 = f2(t(i), c(i), cbar(i));
        k2_2 = f2(t(i) + 0.5*h, c(i) + 0.5*h*k1_2, cbar(i));
        k3_2 = f2(t(i) + 0.5*h, c(i) + 0.5*h*k2_2, cbar(i));
        k4_2 = f2(t(i) + h, c(i) + h*k3_2, cbar(i));
        c(i+1) = c(i) + h/6*(k1_2 + 2*k2_2 + 2*k3_2 + k4_2);
    end

    % Update initial values for next iteration
    c(1) = c(end);
    cbar(1) = cbar(end);

    % Display iteration number and updated initial values
    disp(['Iteration ', num2str(iter), ': c(1) = ', num2str(c(1)), ', cbar(1) = ', num2str(cbar(1))]);
end

Iteration 1: c(1) = 5.8892e-09, cbar(1) = 6.6908
Iteration 2: c(1) = 1.0265e-08, cbar(1) = 4.9741
Iteration 3: c(1) = 1.3516e-08, cbar(1) = 3.6979
Iteration 4: c(1) = 1.593e-08, cbar(1) = 2.7491
Iteration 5: c(1) = 1.7723e-08, cbar(1) = 2.0437
Iteration 6: c(1) = 1.9053e-08, cbar(1) = 1.5194
Iteration 7: c(1) = 2.0039e-08, cbar(1) = 1.1295
Iteration 8: c(1) = 2.0771e-08, cbar(1) = 0.83971
Iteration 9: c(1) = 2.1312e-08, cbar(1) = 0.62426
Iteration 10: c(1) = 2.1712e-08, cbar(1) = 0.46409

plot(t,c)
xlabel('t')
ylabel('c')
set(gca,'Fontsize',16)

Answer 1

[Torsten]

Write a function

function dy = f(t,y)

that accepts t and a vector y of length 2*n with

y = [c_1;...;c_n,cbar_1;...;cbar_n]

and returns

dy = [dc_1/dt;...;dc_n/dt;dcbar_1/dt;...,dcbar_n/dt]

Then write the Runge-Kutta-method of order 4 for a system of differential equations of the form

y' = f(t,y)

and apply it to your above f.

tstart = 0.0;
tend = 5000.0;
h = (tend - tstart)/1000000;
T = (tstart:h:tend).';
n = 10; % Number of stages/bed divisions
Y0 = zeros(1,2*n);
Y0(n+1) = 9;
Y = runge_kutta_RK4(@(t,y)f(t,y,n),T,Y0);
plot(T,Y(:,10))                       % Runge Kutta solution
grid on

function dy = f(t,y,n)

  e = 0.38; % Bed porosity
  rho = 126; % Solvent (CO2) density
  W = 0.95; % CO2 flow rate, kg/h
  V = 0.0004; % Extractor volume, m^3
  Di = 6 * 10^-13; % Diffusion coefficient
  ti = 1736.111; % Internal diffusion time
  c_in = @(t) 0;
  c = y(1:n);
  cbar = y(n+1:2*n);
  dc = zeros(1,n);
  dcbar = zeros(1,n);
  dcbar = -1/ti * (cbar-c/0.2);
  dc(1) = -(W/rho*(c(1)-c_in(t))+(1-e)*V*dcbar(1))/(e*V);
  dc(2:n) = -(W/rho*(c(2:n)-c(1:n-1))+(1-e)*V./(2:n).*dcbar(2:n))./(e*V./(2:n));
  dy = [dc,dcbar];
end

function Y = runge_kutta_RK4(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for i = 2:N
    t = T(i-1);
    y = Y(i-1,:);
    h  = T(i) - T(i-1);
    k0 = f(t,y);
    k1 = f(t+0.5*h,y+k0*0.5*h);
    k2 = f(t+0.5*h,y+k1*0.5*h);
    k3 = f(t+h,y+k2*h);
    Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end

[Jade Pollard]

Your help has been great so far. The curve ideally suppose to replicate this graph shape. Is there an error with our approach ?

[Torsten]

The reason for the discrepancy can only lie in the units (graph seems to be weight-percent), in your inflow setting (c_0(t) which is set to 0), in your vector of initial conditions Y0 and/or in your equations. The numerical method is tested and works.