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
plot(t,c)
xlabel(‘t’)
ylabel(‘c’)
set(gca,’Fontsize’,16)The following is an supercritical CO2 equation of an oil from leaves:

% 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 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
plot(t,c)
xlabel(‘t’)
ylabel(‘c’)
set(gca,’Fontsize’,16) The following is an supercritical CO2 equation of an oil from leaves:

% 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 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
plot(t,c)
xlabel(‘t’)
ylabel(‘c’)
set(gca,’Fontsize’,16) runge kutta, 4th order, supercritical co2 extraction, oil extraction MATLAB Answers — New Questions

Cart

Cart