Ground water contamination due to microbes ( modeling using explicit finite difference method)

The governing equations that I am trying to solve are

subjected to the boundary conditions

Finite difference formulation

Code

% Finite difference grid

dt = 0.1;   %seconds
dz = 0.1;   %cm

T = 3600;  %seconds
L = 70;     %cm

t = linspace(0,T,36000)

% % % 
% % % t = 0:dt:T
% % % z = 0:dz:L

Nt = T/dt;
Nz = L/dz ; %


% Parameters

% Parameters in S

Vw = 0.34 ; % cm/s
D = 1.51e-6 ; % cm2/s
mum = 0.94 ; % microbes/s
ks = 0.9 ; % mg/ml
Y = 6.5 ; % g/mole ATP
rhos = 1 ;
theta = 0.4 ;

%

% Parameters in Ma

ka = 0.85 ;
Mam = 1100 ;
ky = 0.37 ;
b = 4.17e-11 ; % microbes/s

%

% Parameters in M

V = 0.66; % cm/s
tau = 2;

%
kf = 0.95;
Rf = 1 + rhos*kf/theta ;

% Matrix to store solution

S = zeros(Nz, Nt);
Ma = zeros(Nz, Nt);
M = zeros(Nz, Nt);

% Setting initial conditions

S(:, 1) = 10;
Ma(:, 1) = 1000;
M(:, 1) = 1000;

% Setting boundary condition at inlet

S(1, :) = 0;
M(1, :) = 0;





% Solving using explicit method
a1 = (-Vw*dt/(2*dz) + D*dt/(dz*dz))/Rf;
a2 = 1 - 2*D*dt/(Rf*dz*dz);
a3 = (Vw*dt/(2*dz) + D*dt/(dz*dz))/Rf;

c1 = -Vw*dt/(2*dz) - V*dt/(2*dz) + tau*D*dt/(dz*dz);
c2 = 1 - 2*tau*D*dt/(dz*dz);
c3 = Vw*dt/(2*dz) + V*dt/(2*dz) + tau*D*dt/(dz*dz);
for j = 2:Nt
    for i = 2:Nz-1
        % For Solving S

        b1 = mum/((ks + S(i,j-1))*Y);
        b2 = (M(i,j-1) + rhos*Ma(i,j-1)/theta)*dt/Rf;

        S(i,j) = S(i+1,j-1)*a1 + S(i,j-1)*(a2-b1*b2) + S(i-1,j-1)*a3;



        % For Solving M

        d1 = ka*(1 - Ma(i,j-1)/Mam)*dt;
        d2 = (mum*S(i,j-1)/(ks + S(i,j-1)) - b)*dt;
        d3 = (ky*rhos*dt/theta)*Ma(i,j-1);

        M(i,j) = M(i+1,j-1)*c1 + M(i,j-1)*(c2-d1+d2) + M(i-1,j-1)*c3 + d3;
    end

    % Setting Boundary condition at outlet

    i = Nz;

    % For S

    b1 = mum/((ks + S(i,j-1))*Y);
    b2 = (M(i,j-1) + rhos*Ma(i,j-1)/theta)*dt/Rf;

    S(i,j) = S(i-1,j-1)*a1 + S(i,j-1)*(a2-b1*b2) + S(i-1,j-1)*a3;

    % For M

    d1 = ka*(1 - Ma(i,j-1)/Mam)*dt;
    d2 = (mum*S(i,j-1)/(ks + S(i,j-1)) - b)*dt;
    d3 = ky*rhos*dt*Ma(i,j-1);

    M(i,j) = M(i-1,j-1)*c1 + M(i,j-1)*(c2-d1+d2) + M(i-1,j-1)*c3 + d3;

    % Solving Ma

    for i = 1:Nz
        e1 = 1 - ka*theta*M(i,j-1)/(rhos*Mam);
        e2 = mum*S(i,j-1)/(ks + S(i,j-1)) - b;
        e3 = ka*theta*dt*M(i,j-1)/rhos;

        Ma(i,j) = Ma(i,j-1)*(e1 - ky + e2)*dt + e3;
    end
end


%%%%%% Output that i need

plot(t,S(100, :))    

plot(Z,S(:,18000))

plot(t,Ma(100,:))

plot(t,M(100,:))

plot(z,M(:,18000))

The output that is not coming for me but the code is working please tell me the changes that i need to do

Answer 1

[Nivedita]

Hi Haridas!

I understand that you are experiencing an issue while observing the output of your code.

The following changes can help resolve the issues:

• I added this line of code for defining "z" in the given code by referring to your definition of "t" and the comment section below it.

        z = linspace(0,L,700);

As this line of code was not available the code threw an error that the variable "z" was unrecognized.

• To view the different plots in the same figure, I would suggest using "subplot". You can refer to the documentation of "subplot" here: https://www.mathworks.com/help/releases/R2021b/matlab/ref/subplot.html?searchHighlight=subplot&s_tid=doc_srchtitle
• Another observation in the calculations of "M", "S" and "Ma" was that at once instance MATLAB is trying to multiply values which are even greater than the order of e^250, resulting in the answer to be Inf (infinity) as double cannot hold such large values. Since the values calculated in each of the matrices are interdependent on each other and use values from adjacent cells, the next calculation uses Inf resulting in NaN or -Inf.

I hope this answers your question.