Algae growing, concentration curve problem

3 views (last 30 days)
Alfonso on 4 Jan 2024
Edited: Torsten on 6 Jan 2024
Hi, i have made the following code:
function w = Bioreattore_Matlab
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% U(0,t) = 0
% U(L,t) = u (end-1,t)
% U(z, 0) = U0(z)=0
% it is assumed that the concentration of C,N,P are costant
% Problem parameters
alpha1 = 5*10^(-4);
mu0 = 5;
mue = 3.5;
Hl = 5;
T = 8000; % [s]
nz = 100; nt = 2000;
KP = 2 ;
HP= 7;
P = 3 ;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
N = 4;
NC = 2;
KC = 5;
HC = 5 ;
C = 3 ;
PHs3 = 4 ;
PHs4 = 5;
Z = 10;
lambda = 2 ;
rs= 10;
%Check of the problem's parameters
if any([P-PC N-NC] <=0 )
error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T nz-2 nt-2 KP P PC HC KN N NC KC HC C HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz=Z/nz ;
dt = T/nt;
r = alpha1*dt/((dz)^2);
st = 1;
while (2*r) > st
nt = nt+1;
end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
f1= (KP*(P-PC))/(HP+(P-PC));
f2 = (KN*(N-NC))/(HN+(N-NC));
pH = (6.35 - log10(C)) /2 ;
f3 = 1/(1+exp(lambda*(pH-PHs3)));
f4 = 1/(1+exp(lambda*(pH-PHs4)));
%I0 = 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
w(:,1) = 0; %Final concentration of the algae [0]
% Finite-difference method
%for j=2:nt+1
% u(1) = 0; % Condizione al contorno di Dirichlet
% g = (mu0.*I(j))/(Hl+I(j)); decay = g(j)*f1*f2*f3-mue*f4;
% u(2:end-1)=(1-2*r-decay*dt)*u(2:end-1)+r*(u(1:end-2)+u(3:end));
% u(end) = u(end); % Condizione al contorno di Neumann
for j = 1:nt
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j))); %Equation (3) in Algae
%I_in = I0(t(j))*exp(-K*z.').*exp(-rs)*cumtrapz(z.',w(:,j));
g = mu0*I_in./(Hl+I_in); %Equation (2) in Algae
integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j)); %The first part of the equation (1)
%w(1,j+1) = w(2,j+1);
%w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(2:nz-3,j)+w(2:nz-2,j));
w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz,j)+r*(w(1:nz-1,j)+w(3:nz+1,j)); %Equation (1) in Algae
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = 5;%w(nz,j+1); %Concentration at z=0
%w(2:nz,j+1) = w(2:nz,j) + dt*...
%w(2:nz,j+1) = integral_term(j) +
plot(z,flipud(w(:,j+1)),'r*:');
%set(gca,'XDir','reverse');
xlabel('z'); ylabel('w');
title(['t=',num2str(t(j))]);
axis([0 10 0 10]);
pause(0.0001);
end
%plot(z,flipud(w(:,end)))
end
It regards the growing of algae in a tower in which there is light only on 1 side of it. So, first of all, at z=0 we put some seeds so that we have a certain concentration of algae and then when i start the code i would like to see the path of the algae concentration (w) over the whole tower and over the time.
At this point, the code that i've made works really good until T=2000 (seconds) but after the code stops working and so the curve that i see on the printed graph is like "stucked" and so it wont work anymore.
Im loosing my mind to find some problems, but i really dont know what to modify, i've tried to do everything i know. Can you help me to solve si problem, i would like to reach also T=20000 seconds without any problem to my curve printed.
P.S. Attached to this post you find also the pdf file of the article that i read to make the code.
Sam Chak on 5 Jan 2024
Edited: Sam Chak on 5 Jan 2024
I noticed that you used cumtrapz() and trapz() in the code, and I'm checking if you can actually model the differential equations like that. Could you please edit and add comments, such as the equation number, in the code? This will help to verify if the equations can be effectively coded in MATLAB.
Additionally, it seems that the integral terms in Eq. (9) can be approximated by the sum, for i=1 to m=100.
Excerpt: Kumar, K., Pisarenco, M., Rudnaya, M., & Savcenco, V. (2014). A note on analysis and numerics of algae growth. Nonlinear Analysis : Real World Applications, 15(1), 392-403.
Alfonso on 6 Jan 2024
I've modified the code with some comments like Equation (3) in Algae.
Regarding the equation of f1, f2, f3, they are funcion, respectivetely of P, N and C. Now first i want that the code works, so i hypotized P, N and C as constants so alfo f1, f2 and f3 are constant. But while we solve the problem that i explained above, we should go through the implementation of the equation of P, N and C which changes over time and space like you can see in the equations (9) in the algae pdf.

Alan Stevens on 5 Jan 2024
Increase nt to 10000 and it works just fine for T = 20000.
Alfonso on 6 Jan 2024
function w = Bioreattore_Matlab
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% U(0,t) = 0
% U(L,t) = u (end-1,t)
% U(z, 0) = U0(z)=0
% it is assumed that the concentration of C,N,P are costant
% Problem parameters
alpha1 = 5*10^(-4);
mu0 = 5;
mue = 3.5;
Hl = 5;
T = 8000; % [s]
nz = 100; nt = 1000;
KP = 2 ;
HP= 7;
P = 3;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
N = 4;
NC = 2;
KC = 5;
HC = 5 ;
C = 3 ;
PHs3 = 4 ;
PHs4 = 5;
Z = 10;
lambda = 2 ;
rs= 10;
%Check of the problem's parameters
if any([P-PC N-NC] <=0 )
error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T nz-2 nt-2 KP P PC HC KN N NC KC HC C HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz=Z/nz ;
dt = T/nt;
%r = alpha1*dt/((dz)^2);
%st = 1;
%while (2*r) > st
% nt = nt+1;
%end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
%f1= (KP*(P-PC))/(HP+(P-PC));
%f2 = (KN*(N-NC))/(HN+(N-NC));
%pH = (6.35 - log10(C)) /2 ;
%f3 = 1/(1+exp(lambda*(pH-PHs3)));
%f4 = 1/(1+exp(lambda*(pH-PHs4)));
%I0 = 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
w(:,1) = 0; %Final concentration of the algae [0]
% Finite-difference method
%for j=2:nt+1
% u(1) = 0; % Condizione al contorno di Dirichlet
% g = (mu0.*I(j))/(Hl+I(j)); decay = g(j)*f1*f2*f3-mue*f4;
% u(2:end-1)=(1-2*r-decay*dt)*u(2:end-1)+r*(u(1:end-2)+u(3:end));
% u(end) = u(end); % Condizione al contorno di Neumann
for j = 1:nt
f1= (KP*(P(j)-PC))/(HP+(P(j)-PC));
f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
pH = (6.35 - log10(C(j))) /2 ;
f3 = 1/(1+exp(lambda*(pH-PHs3)));
f4 = 1/(1+exp(lambda*(pH-PHs4)));
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
g = mu0*I_in./(Hl+I_in);
integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));
N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j))
P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j))
C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j))
%I_in = I0(t(j))*exp(-K*z.').*exp(-rs)*cumtrapz(z.',w(:,j));
%w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz,j)+r*(w(1:nz-1,j)+w(3:nz+1,j));
w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)
w(1,j+1) = 0%w(2,j+1);
w(nz+1,j+1) = 5;%w(nz,j+1); %Concentration at z=0
%w(2:nz,j+1) = w(2:nz,j) + dt*...
%w(2:nz,j+1) = integral_term(j) +
plot(z,flipud(w(:,j+1)),'b');
%set(gca,'XDir','reverse');
xlabel('z'); ylabel('w');
title(['t=',num2str(t(j))]);
axis([0 10 0 10]);
pause(0.0001);
end
%plot(z,flipud(w(:,end)))
end
Done! But at T=1000 i have problems, after this time you see that the code is instable
Torsten on 6 Jan 2024
Edited: Torsten on 6 Jan 2024
P = 3;
N = 4;
C = 3 ;
Are these really the correct values for P, N and C at time t = 0 or are these different parameters ?
If they are correct, please use them as
P = zeros(nt+1,1)
N = zeros(nt+1,1)
C = zeros(nt+1,1)
P(1) = 3;
N(1) = 4;
C(1) = 3;
before you enter the for-loop.
And as suggested by me and other forum members several times now: choose a smaller time step (thus enlarge nt). Then the stability problems will vanish.
And to avoid the flipping of the solution for w, simply set
w = zeros(nz+1,nt+1);
w(1,:) = 5;
when initializing w and remove the lines
w(1,j+1) = 0%w(2,j+1);
w(nz+1,j+1) = 5;%w(nz,j+1); %Concentration at z=0
within the for-loop.
The condition w(z=Inf) = 0 is not natural - you should take care that the value for w at z = Z never influences what happens within 0 <= z < Z. A condition dw/dz = 0 at z = Z would make more sense.