Algae growing. Concentration curve problem

Question

Alfonso 2024-1-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2073356-algae-growing-concentration-curve-problem

评论： Alfonso 2024-2-2

algae.pdf

Hi! i have made the following code:

function w = Bioreactor
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
alpha1 =  5*10^(-4);                                       % Diffusion coefficient [m^2/s]
mu0 = 5; 
mue = 3.5; 
Hl = 5;
T = 5000;                                                  % [s]
nz = 100; nt = 1000;
KP = 2; 
HP= 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 
lambda = 2 ; 
rs= 10;
P(1)=3;                                                     % Initial chosen concentration of P [g/L]
N(1)=9;                                                     % Initial chosen concentration of N [g/L]
C(1)=30;                                                    % Initial chosen concentration of C [g/L]
Z = 10;                                                     % Lenght of the bioreactor [m]
%Check of the problem's parameters
if  any([P(1)-PC  N(1)-NC] <=0 ) 
    error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T  nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt);                                             % Round nt to an integer
nz = round(nz);                                             % Round nz to an integer
N = zeros(1, nt+1);                                         % Initialize N to zero
P = zeros(1, nt+1);                                         % Initialize P to zero
C = zeros(1, nt+1);                                         % Initialize C to zero
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
% Finite-difference method
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));    % Data organization
    N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j));           % Function of N over the time and space
    P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j));           % Function of P over the time and space
    C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j));           % Function of C over the time and space
    
    w(nz+1,j+1) = 9;                                        % Algae concentration at z=0  [Dirichlet]
    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) = w(2,j+1);                                    % Algae concentration at z=10 [Neumann]
    
    
    plot(z,flipud(w(:,j+1)),'b');
    %set(gca,'XDir','reverse');
    xlabel('z'); ylabel('w');
    title(['t=',num2str(t(j)),' s']); 
    axis([0 10 0 10]);
    pause(0.0001);
end
end

And the code is working perfectly. But the W which is the algae concentration over space and time it has been fixed a value of 9 mg/L at the beginning of this bioreactor, so at z=0, im talking about the following row:

w(nz+1,j+1) = 9;                                        % Algae concentration at z=0  [Dirichlet]

But now i would like to fix this concentration at a certain z which has to be equal to z=Z/2 so at the half of the bioreactor. I know that in order to make the code working for sure i have to modify this row:

    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);

Can you help me to do that? I've tried something but anything really worked.

Thanks, attached to this post there is the file from which i took the equations.

I should have a graph like this:

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Alfonso 2024-1-23

@Torsten

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-1-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2073356-algae-growing-concentration-curve-problem#answer_1395771

编辑：Torsten 2024-1-23

在 MATLAB Online 中打开

Here is a code for a simple heat-conduction equation

dT/dt = d^2T/dz^2
T(z,0) = 0
T(0,t) = 200
T(1,t) = 340
T(0.5,t) = 90

You should be able to adapt it to your application.

L = 1;

zcut = L/2;

n1 = 50;

n2 = 50;

z1 = linspace(0,zcut,n1);

dz1 = z1(2)-z1(1);

z2 = linspace(zcut,L,n2);

dz2 = z2(2)-z2(1);

tstart = 0;

tend = 1;

dt = 1e-5;

nt = round((tend-tstart)/dt);

t = linspace(tstart,tend,nt);

T_at_zcut = 90;

T_at_0 = 200;

T_at_L = 340;

T = zeros(n1+n2-1,nt);

T(1,:) = T_at_0;

T(end,:) = T_at_L;

T(n1,:) = T_at_zcut;

for i = 1:nt-1

T(2:n1-1,i+1) = T(2:n1-1,i) + dt/dz1^2*(T(3:n1,i)-2*T(2:n1-1,i)+T(1:n1-2,i));

T(n1+1:n1+n2-2,i+1) = T(n1+1:n1+n2-2,i) + dt/dz2^2*(T(n1+2:n1+n2-1,i)-...

2*T(n1+1:n1+n2-2,i)+T(n1:n1+n2-3,i));

end

plot([z1,z2(2:end)],[T(:,1),T(:,100),T(:,250),T(:,500),T(:,750),T(:,1000),T(:,2500),T(:,5000),T(:,7500),T(:,end)])

26 个评论
显示 24更早的评论隐藏 24更早的评论

Alfonso 2024-1-24

在 MATLAB Online 中打开

Ohh yes, i got this code. What do you think? I think it's working!

function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
alpha1 = 0.005;
mu0 = 5; 
mue = 3.5; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Check of the problem's parameters
if  any([P(1)-PC  N(1)-NC] <= 0 ) 
    error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
% Initialization
z = linspace(0, Z, nz+1);
t = linspace(0, T, nt+1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt+1);
P = zeros(1, nt+1);
C = zeros(1, nt+1);
I0 = @(t) 100*max(sin((t-0.4*3600)/(6.28)), sin(0));
w = zeros(nz+1, nt+1);
w(50,:) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1+n2-1, nt+1);
% Finite-difference method
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));
    
    w(nz+1,j+1) = w(nz,j+1);
    w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + alpha1 * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
    w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + alpha1 * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);
    
    plot(z,w(:,j+1),'b');
    xlabel('z'); ylabel('w');
    title(['t=',num2str(t(j)),' s']); 
    axis([0 10 0 10]);
    pause(0.0001);
end
end

Do you think that the results can be associated to the curve that i draw? I think so

Alfonso 2024-1-25

编辑：Alfonso 2024-1-26

在 MATLAB Online 中打开

Ok @Torsten, i modified the code as you told me and so i get this code:

function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5; 
mue = 2; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Check of the problem's parameters
if any([P(1) - PC  N(1) - NC] <= 0 ) 
    error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);
% Finite-difference method
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));
    
    w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
    w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = w(nz,j+1);
    
    plot(z, w(:, j + 1), 'b');
    xlabel('z'); ylabel('w');
    title(['t=', num2str(t(j)), ' s']); 
    axis([0 10 0 10]);
    pause(0.0001);
end
end

As a plot i have for example:

That seems correct because from the center value, i see 2 different curves with 2 different values on the borders. But, the I_in, which is the irradiation, is it put only on 1 side? Because at least i see that the 2 curves are increasing a lot, yes with 2 different values of w, but they are increasing together and i don't want something like this. In this case the plot should have a curve that has a permanent lower values of W and the other curve which is the oppost, i mean this other curve has to increase the w value with a different velocity due to the fact that the algae are irradiated by the light.

Alfonso 2024-1-31

在 MATLAB Online 中打开

@Torsten

I have the following updated code:

function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5; 
mue = 2; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Check of the problem's parameters
if any([P(1) - PC  N(1) - NC] <= 0 ) 
    error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);
% Finite-difference method
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));
    
    w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
    w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = w(nz,j+1);
    
    plot(z, w(:, j + 1), 'b');
    xlabel('z'); ylabel('w');
    title(['t=', num2str(t(j)), ' s']); 
    axis([0 10 0 10]);
    pause(0.0001);
end
end

I understood that I_in is fixed for every value of the z axis, but instead i would like to have it only for like z=10 for every J value. I've tried to modify it and i got:

function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5; 
mue = 2; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Check of the problem's parameters
if any([P(1) - PC  N(1) - NC] <= 0 ) 
    error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);
I_in = zeros(nz + 1, nt + 1);
    for j = 1:nt
        I_in(:, j) = I0(t(j)) * exp(-K * z.') .* exp((-rs) * cumtrapz(z, w(:, j).')');
    end
% Finite-difference method
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_10 = I_in(10, j);
    g = mu0 * I_in_10 / (Hl + I_in_10);
    
    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));
    
    w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
    w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = w(nz,j+1);
    
    plot(z, w(:, j + 1), 'b');
    xlabel('z'); ylabel('w');
    title(['t=', num2str(t(j)), ' s']); 
    axis([0 10 0 10]);
    pause(0.0001);
end
end

The problem is that the code works but doesn't give me anything as result. And i dont know why D:

Can you help me?

Torsten 2024-1-31

在 MATLAB Online 中打开

Bioreactor_central()

function Bioreactor_central

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters

Dm = 5 * 10^(-4);

mu0 = 5;

mue = 2;

Hl = 5;

T = 5000;

nz = 100; %nt = 1000;

KP = 2;

HP = 7;

PC = 1;

K = 2;

KN = 3;

HN = 14.5 * 10^(-6);

NC = 2;

KC = 5;

HC = 5 ;

PHs3 = 4 ;

PHs4 = 5;

lambda = 2 ;

rs = 10;

P(1) = 3;

N(1) = 9;

C(1) = 30;

Z = 10;

% Stability

dz = Z / nz;

st = 1 / 2;

dt = st * dz^2 / Dm;

nt = T / dt;

% Check of the problem's parameters

if any([P(1) - PC N(1) - NC] <= 0 )

error('Check problem parameters')

end

if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)

error('Check problem parameters')

end

% Initialization

z = linspace(0, Z, nz + 1);

t = linspace(0, T, nt + 1);

nt = round(nt);

nz = round(nz);

N = zeros(1, nt + 1);

P = zeros(1, nt + 1);

C = zeros(1, nt + 1);

I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));

w = zeros(nz + 1, nt + 1);

w(50, :) = 5;

I_in = zeros(nz + 1,nt);

% Finite-difference method

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(:,j) = I0(t(j)) .* exp(-K * z.') .* exp((-rs) .* cumtrapz(z.', w(:, j)));

I_in_10 = I_in(Z,j);

g_in_10 = mu0 * I_in_10 / (Hl + I_in_10);

integral_term = f1 * f2 * f3 * trapz(z, g_in_10 .* 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));

w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g_in_10 + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);

w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g_in_10 + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);

w(1,j+1) = w(2,j+1);

w(nz+1,j+1) = w(nz,j+1);

end

plot(z,w(:,end))

end

Alfonso 2024-2-1

在 MATLAB Online 中打开

@Torsten

I've tried to split everything and putting I_in equal to 0 from z=0 to z=5, so that i have the following code:

function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5; 
mue = 2; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
P_1(1) = 3;
N_1(1) = 9;
C_1(1) = 30;
Z = 10;
% Check of the problem's parameters
if any([P(1) - PC  N(1) - NC] <= 0 ) 
    error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
N_1 = zeros(1, nt + 1);
P_1 = zeros(1, nt + 1);
C_1 = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);
% Finite-difference method
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))); 
    f11 = (KP * (P_1(j) - PC)) / (HP + (P_1(j) - PC));
    f22 = (KN * (N_1(j) - NC)) / (HN + (N_1(j) - NC));
    pH1 = (6.35 - log10(C_1(j))) / 2 ;
    f33 = 1 / (1 + exp(lambda * (pH1 - PHs3)));
    f44 = 1 / (1 + exp(lambda * (pH1 - PHs4))); 
    
    I_in = zeros(nz + 1, 1);
    I_in_1_5 = 0;%I0(t(j)) .* exp(-K * z(1:5).') .* exp((-rs) .* cumtrapz(z(1:5).', w(1:5, j)));
    I_in_5_10 = I0(t(j)) .* exp(-K * z(5:10).') .* exp((-rs) .* cumtrapz(z(5:10).', w(5:10, j)));
    g_1_5 = mu0 * I_in_1_5 ./(Hl + I_in_1_5);
    g_5_10 = mu0 * I_in_5_10 ./(Hl + I_in_5_10);
    integral_term_1_5 = f11 * f22 * f33 * trapz(z(1:5), g_1_5 .* w(1:5, j));
    integral_term_5_10 = f1 * f2 * f3 * trapz(z(5:10), g_5_10 .* w(5:10, j));
    
    N_1(j + 1) = N(j) + dt * (-1 / Z * integral_term_1_5 * N(j));
    P_1(j + 1) = P(j) + dt * (-1 / Z * integral_term_1_5 * P(j));
    C_1(j + 1) = C(j) + dt * (-1 / Z * integral_term_1_5 * C(j));
    
    N(j + 1) = N(j) + dt * (-1 / Z * integral_term_5_10 * N(j));
    P(j + 1) = P(j) + dt * (-1 / Z * integral_term_5_10 * P(j));
    C(j + 1) = C(j) + dt * (-1 / Z * integral_term_5_10 * C(j));
    
    w(2:49,j+1) = w(2:49,j) + dt * (f11*f22*f33*w(2:49,j).*g_1_5(2:49) + Dm * (w(3:50,j) - 2*w(2:49,j) + w(1:48,j))/dz^2) - mue*f4*w(2:49,j);
    w(51:nz,j+1) = w(51:nz,j) + dt * (f1*f2*f3*w(51:nz,j).*g_5_10(51:nz) + Dm * (w(52:nz+1,j)-2*w(51:nz,j)+w(50:nz-1,j))/dz^2) - mue*f4*w(51:nz,j);
    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = w(nz,j+1);
    
    plot(z, w(:, j + 1), 'b');
    xlabel('z'); ylabel('w');
    title(['t=', num2str(t(j)), ' s']); 
    axis([0 10 0 10]);
    pause(0.0001);
end
%plot(z,w(end));
end

But i have the following error:

Index exceeds the number of array elements. Index must not exceed 1.
Error in Bioreactor_central1 (line 110)
    w(2:49,j+1) = w(2:49,j) + dt * (f11*f22*f33*w(2:49,j).*g_1_5(2:49) + Dm * (w(3:50,j) - 2*w(2:49,j) + w(1:48,j))/dz^2) - mue*f4*w(2:49,j);

Do you know how can i solve?

Alfonso 2024-2-1

在 MATLAB Online 中打开

You are perfectly right! I've returned to the old code, the following:

function w = Bioreattore
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
alpha1 =  5*10^(-4); 
mu0 = 5; 
mue = 3.5; 
Hl = 5;
T = 10000;                                                  % [s]
nz = 100; nt = 1000;
KP = 2 ; 
HP= 7; 
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 
lambda = 2 ; 
rs= 10;
P(1)=3;                                                     % Initial chosen concentration of P
N(1)=9;                                                     % Initial chosen concentration of N 
C(1)=30;                                                    % Initial chosen concentration of C
Z = 10;                                                     % Lenght of the bioreactor
%Check of the problem's parameters
if  any([P(1)-PC  N(1)-NC] <=0 ) 
    error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T  nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt);                                             % Round nt to an integer
nz = round(nz);                                             % Round nz to an integer
N = zeros(1, nt+1);                                         % Initialize N to zero
P = zeros(1, nt+1);                                         % Initialize P to zero
C = zeros(1, nt+1);                                         % Initialize C to zero
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
% Finite-difference method
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));    % Data organization
    N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j));           % Function of N over the time and space
    P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j));           % Function of P over the time and space
    C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j));           % Function of C over the time and space
    
    w(nz+1,j+1) = 5;                                        % Algae concentration at z=0    
    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)- mue*f4*w(2:nz,j);
    w(1,j+1) = w(2,j+1);                                    % Algae concentration at z=10
   
    
    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
end

But at least while i remove put I_in equal to 0 i have the same plot, as follows:

What it's you advise?

Alfonso 2024-2-2

在 MATLAB Online 中打开

@Torsten

Ok, so. I did it, the code is the following:

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 = 10000; % [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;
P(1)=3;
N(1)=9;
C(1)=30;
%Check of the problem's parameters
if  any([P(1)-PC  N(1)-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/2;
dt = st*dz^2/alpha1;
nt = T/dt;
%while (2*r) > st
 %   nt = nt+1;
%end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt); % Arrotonda nt a un numero intero
N = zeros(1, nt+1); % Inizializzazione di N a zero
P = zeros(1, nt+1); % Inizializzazione di P a zero
C = zeros(1, nt+1); % Inizializzazione di C a zero
%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));
nz = round(nz);
nt = round(nt);
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) = 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(I_in),'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

But the plots are not good as i expected, look at that:

They are the only plot that i see for all the time t is like 1 second the first plot and the next second the other. What do you think? What's you advise, they are just the I_in plot, dont read on yaxes the label.

Alfonso 2024-2-2

在 MATLAB Online 中打开

@Torsten

Regarding the old code:

function w = Bioreattore
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
alpha1 =  5*10^(-4); 
mu0 = 0; 
mue = 99999999;%3.5; 
Hl = 5;
T = 15000;                                                  % [s]
nz = 100; nt = 1000;
KP = 2 ; 
HP= 7; 
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 
lambda = 2 ; 
rs= 10;
P(1)=3;                                                     % Initial chosen concentration of P
N(1)=9;                                                     % Initial chosen concentration of N 
C(1)=30;                                                    % Initial chosen concentration of C
Z = 10;                                                     % Lenght of the bioreactor
%Check of the problem's parameters
%if  any([P(1)-PC  N(1)-NC] <=0 ) 
%    error('Check problem parameters')
%end
%if any([alpha1 mu0 Hl Z T  nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
%    error('Check problem parameters')
%end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt);                                             % Round nt to an integer
nz = round(nz);                                             % Round nz to an integer
N = zeros(1, nt+1);                                         % Initialize N to zero
P = zeros(1, nt+1);                                         % Initialize P to zero
C = zeros(1, nt+1);                                         % Initialize C to zero
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
% Finite-difference method
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));     % Data organization
    N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j));           % Function of N over the time and space
    P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j));           % Function of P over the time and space
    C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j));           % Function of C over the time and space
    
    w(nz+1,j+1) = 5;                                         % Algae concentration at z=0    
    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)- mue*f4*w(2:nz,j);
    w(1,j+1) = w(2,j+1);                                     % Algae concentration at z=10
   
    
    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
end

I see that if i change every parameter it doesn't change anything in the curve, itìs not normal, it seems that the w equation doesnt read the initial parameter, it gives me everytime the same curve also if i put everything 0. I just need a code that reads those changes at least, can you check?

Alfonso 2024-2-2

Ok, man i will check. At least i didn't check the dimension before D:

Alfonso 2024-2-5

在 MATLAB Online 中打开

@Torsten also with the right dimension the code doesn't work very well. Can i ask you a question? If i have this code:

function w = Bioreattore
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm =  5*10^(-4);                                            %m^2/s              
mu0 = 8000%2.5;                                                  % L/day
mue = 3.5;                                                  % 1/day
Hl = 5;     %half-saturation intensity W/(m^2*day)
T = 10000;                                                  % [s]
nz = 100; nt = 1000;
KP = 2 ; 
HP= 7;             %half-saturation concentrations of Phosphates
PC = 1;     %critical concentration of the phosphates dopo la quale la crescita di w va a 0.
K = 2;
KN = 3 ;
HN = 0.02;   %14.5*10^(-6);  %half-saturation concentrations of nitrates [mgN/L]
NC = 2;     %critical concentration of the nitrates dopo la quale la crescita di w va a 0.
KC = 5;
HC = 5 ;
PHs3 = 9 ; %describes the "switching" value of pH at which the growth increases if all other factors are kept unchanged.
 
lambda = 2 ;                %sharpness of the profile of the f3 function
rs= 10;                                                     %L*m/d
P(1)=50;                                                     % Initial chosen concentration of P [mg/L]
N(1)=330;                                                     % Initial chosen concentration of N [mg/L]
C(1)=2100;                                                    % Initial chosen concentration of C/ Aumentando C, si abbassa il pH(+acido) e aumenta la morte algale [mg/L]
Z = 10;                                                     % Lenght of the bioreactor
%Check of the problem's parameters
if  any([P(1)-PC  N(1)-NC] <=0 ) 
    error('Check problem parameters')
end
if any([Dm mu0 Hl Z T  nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/Dm;
nt = T/dt;
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt);                                             % Round nt to an integer
nz = round(nz);                                             % Round nz to an integer
N = zeros(1, nt+1);                                         % Initialize N to zero
P = zeros(1, nt+1);                                         % Initialize P to zero
C = zeros(1, nt+1);                                         % Initialize C to zero
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
% Finite-difference method
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   %relation between the pH and the concentration of carbon dioxide
    f3 = 1/(1+exp(lambda*(pH-PHs3)))   %The growth rate dependence function. Si riduce se il tutto diventa più acido (pH scende)..
    
    Ha = mue*f3;       %death rate of the algae. [gCOD/L]
    
    I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));     % W/m^2*s
    g = mu0*I_in./(Hl+I_in);
    
    
    integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));     % Data organization
    N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j));           % Function of N over the time and space
    P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j));           % Function of P over the time and space
    C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j));           % Function of C over the time and space
    
    w(nz+1,j+1) = 5;                                         % Algae concentration at z=0    [gSST/L]
    w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + Dm * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)- Ha*w(2:nz,j);
    w(1,j+1) = w(2,j+1);                                     % Algae concentration at z=10
   
    
    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
end

How can i modify the plot to obtain the following graphs:

请先登录，再进行评论。

Algae growing. Concentration curve problem

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

26 个评论
显示 24更早的评论隐藏 24更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Algae growing. Concentration curve problem

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

26 个评论 显示 24更早的评论隐藏 24更早的评论

另请参阅

类别

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

26 个评论
显示 24更早的评论隐藏 24更早的评论