Solving discretized algebraic equations using nested FOR loops

Question

Zain Shabeeb 2019-7-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/471919-solving-discretized-algebraic-equations-using-nested-for-loops

编辑： Zain Shabeeb 2019-7-16

Hi, I am trying to model a process on MATLAB which has algebraic equations and is also time dependent. I want to discretize the d/dt term using a 'for' loop, making the whole model only algebraic. The reason I am not using an 'ode' solver is because I want to optimize this model on another software called GAMS which does not deal with differential equations.

I am having to use nested 'for' loops because the variables depend on 2 directions ('m' direction and 'n' direction), and time, therefore, 3 nested 'for' loops. I have attached the model equations.

I believe one problem is with my nested for loops. There is no error, however the 'for' loop does not assign values to each cell like it should (i.e. the for loop does not do anything).

Furthermore, I believe there is also a problem with the discretization. I have used the backward differential method as suggested by the authors of the model. the DqDt discretization gives 'NaN' for all cells from t=2:tf.

I would be grateful for any help regarding this problem. The model seems lengthy, but I believe that the problem is limited to the nested 'for' loops, because everything before them seems to work okay.

Thanks in advance.

clear all
%% Parameters
    % Molar conductivities (S.m^2.mol^-1)
        % Molar conductivities in liquid phase:
LambdaH_L = 0.035;                    % Molar conductivity of H+ in liquid phase
LambdaNa_L = 0.005;                   % Molar conductivity of Na+ in liquid phase
LambdaA_L = 0.004;                    % Molar conductivity of A- in liquid phase
        % Molar conductivities in solid phase:
LambdaH_S = 0.00192;                  % Molar conductivity of H+ in solid phase
LambdaNa_S = 0.00059;                 % Molar conductivity of Na+ in solid phase
    % Mobilities (m^2.V^-1.s^-1)
        % Mobilities in liquid phase:
uH_L = 36*10^-8;                      % Mobility of H+ in liquid phase
uNa_L = 5.19*10^-8;                   % Mobility of Na+ in liquid phase
uA_L = 4.15*10^-8;                    % Mobility of A- in liquid phase
        
        % Mobilities in solid phase:
uH_S = 19.9*10^-9;                    % Mobility of H+ in solid phase
uNa_S = 6.1*10^-9;                    % Mobility of Na+ in solid phase
    % General system parameters:
Q = 1790;                             % Capacity of ion exchange resin in eq/m^3
K = 1.56;                             % Ion exchange equilibrium constant for Na+/H+
epsilon = 0.4;                        % Bed porosity
Cfeed = 100;                          % Concentration of feed (mol/m^3)
L = 0.015;                            % Bed length = 15mm
W = 0.015;                            % Bed width = 15mm
S = 0.002;                            % Contact area of each membrane = 20cm^2
F = 3.27*10^-7;                       % Feed flowrate = 19.6mL/min
zNa = 1;                              % Valence of Na+ ion
zH = 1;                               % Valence of H+ ion
zA = -1;                              % Valence of A- ion
zA_mod = 1;                           % Positive valence of A- ions                        
a = 0.63;                             % Parameter a for Kcell equation  
b = 0.01;                             % Parameter b for Kcell equation  
c = 0.34;                             % Parameter c for Kcell equation  
d = 0.95;                             % Parameter d for Kcell equation  
e = 0.05;                             % Parameter e for Kcell equation  
V = 0.0001;                           % Volume of bed in m^3
U = 2;                                % Constant Voltage in Volts
%% Defining steps
tf = 10000;                           % Total time = 100s
dt = 100;                             % One time step = 0.1s
n = [0:0.001:L];                    % Vector of  steps in vertical direction
m = [0:0.001:W];                    % Vector of steps in horizontal direction
t = [0:dt:tf];                      % Vector of time steps
N = numel(n);                       % Number of steps in the vertical direction
M = numel(m);                       % Number of steps in the horizontal direction
T = numel(t);                       % Number of time steps
%% Parameters depending only on N and M
Sc = S/N;                           % Contact area of a cell
Fc = F/M;                           % Flowrate in each cell
Wc = W/M;                           % Cell width
Lc = L/N;                           % Cell length
Vc = V/(N*M);                       % Cell volume
%% Initializing terms before for loop
    % Initializing variables
CNa = zeros(N,M,T);                 % Concentration of Na+ in the liquid phase initialized before the for loop
CH = zeros(N,M,T);                  % Concentration of H+ in the liquid phase initialized before the for loop
CA = zeros(N,M,T);                  % Concentration of A+ in the liquid phase initialized before the for loop
CHA = zeros(N,M,T);                 % Concentration of HA in the liquid phase initialized before the for loop
CAtot = zeros(N,M,T);
qNa = zeros(N,M,T);                 % Concentration of Na+ in the solid phase initialized before the for loop
qH = zeros(N,M,T);                  % Concentration of H+ in the solid phase initialized before the for loop
Kcell_L = zeros(N,M,T);             % Cell conductivity in the liquid phase initialized before the for loop
Kcell_S = zeros (N,M,T);            % Cell conductivity in the solid phase initialized before the for loop
Kcell = zeros(N,M,T);               % Cell conductivity initialized before the for loop
tNa = zeros(N,M,T);                 % Transport number of Na+ initialized before the for loop
tA = zeros(N,M,T);                  % Transport number of A- initialized before the for loop
tH = zeros(N,M,T);                  % Transport number of H+ initialized before the for loop
Cond_cell = zeros(N,M,T);           % Cell conductance initialized before the for loop
Resistance_cell = zeros(N,M,T);     % Cell resistance initialized before the for loop
Iline_inv = zeros(N,T);             % 1/Iline(n) initialized before the for loop
Iline = zeros(N,T);                 % Iline(n) initialized before the for loop
Iline_mult = zeros(N,M,T);          % Iline array for multiplication in mass balance equations (cannot use Iline for this because of array dimensions)
   
     %Initializing time differential terms
DCDt_Na = zeros(N,M,T);             % DCDt in Na+ mass balance first accumulation  term                 
DCDt_Atot = zeros(N,M,T);           % DCDt in Total acid mass balance accumulation term
DCDt_H = zeros(N,M,T);              % DCDt in H+ mass balance first accumulation  term 
DqDt_Na = zeros(N,M,T);             % DqDt in Na+ mass balance second accumulation term
DqDt_H = zeros(N,M,T);              % DqDt in H+ mass balance second accumulation term
      
%% Initial conditions at t = 0 (Na form)
    % Initial conditions for variables
CNa(:,:,1) = Cfeed;                 % Concentration of Na+ ions in liquid phase at t = 0 
CH(:,:,1) = 0;                      % Concentration of H+ ions in liquid phase at t = 0
CA(:,:,1) = CNa(:,:,1);             % Concentration of A- ions in liquid phase at t = 0 
CHA(:,:,1) = 0;                     % Concentration of HA in the liquid phase at t = 0
CAtot(:,:,1) = CNa(:,:,1);         % Concentration of total acid in liquid phase at t = 0
qNa(:,:,1) = Q;                     % Concentration of Na+ ions in solid phase at t = 0
qH(:,:,1) = 0;                      % Concentration of H+ ions in solid phase at t = 0
Kcell_L(:,:,1) = zNa*CNa(:,:,1)*LambdaNa_L + zH*CH(:,:,1)*LambdaH_L + zA_mod*CA(:,:,1)*LambdaA_L;                                   % Cell conductivity in the liquid phase at t = 0
Kcell_S(:,:,1) = zNa*qNa(:,:,1)*LambdaNa_S + zH*qH(:,:,1)*LambdaH_S;                                                                % Cell conductivity in the solid phase at t = 0
Kcell(:,:,1) = ((a*Kcell_L(:,:,1).*Kcell_S(:,:,1))./(d*Kcell_L(:,:,1) + e*Kcell_S(:,:,1))) + b*Kcell_S(:,:,1) + c*Kcell_L(:,:,1);     % Cell conductivity at t = 0
tNa(:,:,1) = (zNa*uNa_S*qNa(:,:,1) + zNa*uNa_L*CNa(:,:,1))./(zNa*uNa_S*qNa(:,:,1) + zNa*uNa_L*CNa(:,:,1) + zH*uH_S*(Q-qNa(:,:,1)) + zH*uH_L*CH(:,:,1) + zA_mod*uA_L*CA(:,:,1)); % Transport number of Na+ at t = 0
tA(:,:,1) = (zA_mod*uA_L*CA(:,:,1))./(zNa*uNa_S*qNa(:,:,1) + zNa*uNa_L*CNa(:,:,1) + zH*uH_S*(Q-qNa(:,:,1)) + zH*uH_L*CH(:,:,1) + zA_mod*uA_L*CA(:,:,1));                        % Transport number of A- at t = 0
tH(:,:,1) = (zH*uH_S*qH(:,:,1) + zH*uH_L*CH(:,:,1))./(zNa*uNa_S*qNa(:,:,1) + zNa*uNa_L*CNa(:,:,1) + zH*uH_S*(Q-qNa(:,:,1)) + zH*uH_L*CH(:,:,1) + zA_mod*uA_L*CA(:,:,1));        % Transport number of A- at t = 0
Cond_cell(:,:,1) = Kcell(:,:,1)*(Sc/Wc);            % Cell conductance at t = 0
Resistance_cell(:,:,1) = 1./(Cond_cell(:,:,1));     % Cell resistance at t = 0
Iline_inv(:,1) = sum(Resistance_cell(:,:,1),2);     % Inverse of current of line at t = 0
Iline(:,1) = 1./Iline_inv(:,1);                     % Current of a line at t = 0
for i = 1:M
    Iline_mult(i,:,1) = Iline(i,1);
end
 
CAtot(:,:,1) = CA(:,:,1) + CHA(:,:,1);
CA(:,:,1) = CH(:,:,1) + CNa(:,:,1);
qNa(:,:,1) = (Q*K*CNa(:,:,1))./(CH(:,:,1) + K*CNa(:,:,1));
qH(:,:,1) = Q - qNa(:,:,1);
%% Boundary Conditions
    % Boundary conditions 1, for m = 1, n = 2:N
for i = 2:N 
    CNa(i-1,1,:) = (Fc*CNa(i,1,:) + (tNa(i,1,:).*Iline_mult(i,1,:))/(96485.3*zNa) + Vc*epsilon*DCDt_Na(i,1,:) + Vc*DqDt_Na(i,1,:))./Fc;
    CH(i-1,1,:) = (((-Iline_mult(i,1,:))/(96485.3*zH)) + Fc*CH(i,1,:) + (tH(i,1,:).*Iline_mult(i,1,:))./(96485.3*zH) + Vc*epsilon*DCDt_H(i,1,:) + Vc*DqDt_H(i,1,:))./Fc; 
    CAtot(i-1,1,:) = (((-tA(i,2,:).*Iline_mult(i,2,:))/(96485.3*zA)) + Fc*CA(i,1,:) + Vc*epsilon*DCDt_Atot(i,1,:))./Fc;
    
    CAtot(i-1,M,:) = (Fc*CA(i,M,:) + ((tA(i,M,:).*Iline_mult(i,M,:))/(96485.3*zA))+ Vc*epsilon*DCDt_Atot(i,M,:))./Fc;
    
end
 % Boundary conditions 2, for n = 1, m = 1:M
 
for j = 1:M 
     CA(1,j,:) = CNa(1,j,:);
     CH(1,j,:) = 0;
end
   
   
%% For loops to assign values to each cell
 
for i = 2:N 
    for j = 2:M-1 
        for k = 2:dt:T
            
            DCDt_Na(i,j,k) = (CNa(i,j,k)-CNa(i,j,k-1))./(t(k)-t(k-1));
            DqDt_Na(i,j,k) = (qNa(i,j,k)-qNa(i,j,k-1))./(t(k)-t(k-1));
            DCDt_Atot(i,j,k) = (CAtot(i,j,k)-CAtot(i,j,k-1))./(t(k)-t(k-1));
            DCDt_H(i,j,k) = (CH(i,j,k)-CH(i,j,k-1))./(t(k)-t(k-1));
            DqDt_H(i,j,k) = (qH(i,j,k)-qH(i,j,k-1))./(t(k)-t(k-1));
    
            CNa(i-1,j,k) = (((-tNa(i,j-1,k).*Iline_mult(i,j-1,k))/(96485.3*zNa)) + Fc*CNa(i,j,k) + ((tNa(i,j,k).*Iline_mult(i,j,k))/(96485.3*zNa)) + Vc*epsilon*DCDt_Na(i,j,k) + Vc*DqDt_Na(i,j,k))./Fc;
            CAtot(i-1,j,k) = (((-tA(i,j+1,k).*Iline_mult(i,j+1,k))/(96485.3*zA)) + Fc*CA(i,j,k) + ((tA(i,j,k).*Iline_mult(i,j,k))/(96485.3*zA)) + Vc*epsilon*DCDt_Atot(i,j,k))./Fc;
            CHA(i,j,k) = CAtot(i,j,k) - CA(i,j,k);
            CA(i,j,k) = CH(i,j,k) + CNa(i,j,k);
            qNa(i,j,k) = (Q*K*CNa(i,j,k))./(CH(i,j,k) + K*CNa(i,j,k));
            qH(i,j,k) = Q - qNa(i,j,k);
            Kcell_L(i,j,k) = zNa*CNa(i,j,k)*LambdaNa_L + zH*CH(i,j,k)*LambdaH_L + zA_mod*CA(i,j,k)*LambdaA_L;                                   
            Kcell_S(i,j,k) = zNa*qNa(i,j,k)*LambdaNa_S + zH*qH(i,j,k)*LambdaH_S;                                                               
            Kcell(i,j,k) = ((a*Kcell_L(i,j,k).*Kcell_S(i,j,k))./(d*Kcell_L(i,j,k) + e*Kcell_S(i,j,k))) + b*Kcell_S(i,j,k) + c*Kcell_L(i,j,k);     
            tNa(i,j,k) = (zNa*uNa_S*qNa(i,j,k) + zNa*uNa_L*CNa(i,j,k))./(zNa*uNa_S*qNa(i,j,k) + zNa*uNa_L*CNa(i,j,k) + zH*uH_S*(Q-qNa(i,j,k)) + zH*uH_L*CH(i,j,k) + zA_mod*uA_L*CA(i,j,k)); 
            tA(i,j,k) = (zA_mod*uA_L*CA(i,j,k))./(zNa*uNa_S*qNa(i,j,k) + zNa*uNa_L*CNa(i,j,k) + zH*uH_S*(Q-qNa(i,j,k)) + zH*uH_L*CH(i,j,k) + zA_mod*uA_L*CA(i,j,k));                       
            tH(i,j,k) = (zH*uH_S*qH(i,j,k) + zH*uH_L*CH(i,j,k))./(zNa*uNa_S*qNa(i,j,k) + zNa*uNa_L*CNa(i,j,k) + zH*uH_S*(Q-qNa(i,j,k)) + zH*uH_L*CH(i,j,k) + zA_mod*uA_L*CA(i,j,k));       
            Cond_cell(i,j,k) = Kcell(i,j,k)*(Sc/Wc);            
            Resistance_cell(i,j,k) = 1./(Cond_cell(i,j,k));    
            Iline_inv(i,k) = sum(Resistance_cell(i,j,k),2);    
            Iline(i,k) = 1./Iline_inv(i,k);  
            for l = 1:M
                 Iline_mult(l,:,k) = Iline(l,k);
            end
            
        end
    end
end