Unknown For Loop Placement Error / Increment Error

Hello, I have a script for modeling a CO2 Scrubbing system found in smoke stacks of power plants which computes the concentration of CO2 using the finite difference method. The script iteratively computes the concentration at a (z,r) coordinate and compares it to a local index of (i,j). I am being asked to compute the concentration of the exhaust gas for a range of smoke stack heights. The script successfully solves for the concentration when individual values of length are used, however, I would like to solve for 50 different concentrations with 50 different lengths. To do this, I took the section of code which workes for individual values, and encased it in a for loop which updates the length after the concentration gradient has been determined. My problem is when running the script, it appears that only the first value of the length is used, yet it still returns 50 values of concentrations (they are all the same since the first length is being utilized). I would like to know if anyone can recognize where the issuse might occur, potentially from loop placement / placement of where the length value is incrimented.

***Additionally, the script is specifically looking to calculate the concentration at the exit of the smoke stack which is represented as C_TOP. This does not cause addtional problems with computation, but can be used with C_avg to see that the results are only evaluate for L(ii)=L(1).

Thanks for any help,

Kyle

%% System inputs %%

%Operational Parameters:
C0 = 0.3;        %Inlet mass fraction [-]
Cw = 0.0;        %Reference mass fraction with amine [-]
kw = 1;          %Convective mass transfer coefficent [m/s]
Q  = 1;          %Volumetric flow rate [m3/s]
D  = 0.02;       %Diffusion coefficent (CO2 - amine) [m2/s]


%Design Parameters:
L = linspace(0.1,20,50);        %length [m] between 0 and 20 meters
R = 0.8;         %radius [m] between 0.8 and 1.6 meters
ax = Q/(pi*R^2); %Representative advective velocity

%Design Assessment Parameters:
C_target = 0.1;     %target concentration
Vs0      = 20;      %cost of scrubber + operational costs
ms       = 1.4;
mc       = 100;



%% Model Solution %% 

% numerical:
N = 50;        %number of nodes
Nr = N;        %number of nodes along r
Nz = N;        %number of nodes along z

% tolerance tolerance for solution
tolerance = 1.0e-4;
% define a 2D grid with nodal values of variable:
C = zeros( Nz, Nr, N );

% set up: ---------------------------------------------------------------
for ii = 1: N       %start of for loop used to incrment L
    
%everything from here down works for individual values of L and R but in
%that case L = L, not L = L(ii)

% spacing between nodes in each direction:
deltar = R / ( Nr - 1 );
deltaz = L(ii) / ( Nz - 1 );       %this is where the length value is update for each itteration of for loop (or so i think)

% solution procedure: ---------------------------------------------------
Cold(:,:,ii)    = C(:,:,ii);       % auxiliary array to store solution
error   = 1.0e+8;                  % initialize error as a large number
nite    = 0;                       % iteration counter
nitemax = 1e+5;                    % maximum number of iterations

    % iterate until convergence:
    while ( error > tolerance ) && ( nite < nitemax )

    % update iteration counter:
    nite = nite + 1;

    % estimate new value in nodes from the old values:
        for i = 1 : Nz
            for j = 1 : Nr
            
            % define coordinates for point with indexes (i, j):
            r = deltar * ( j - 1 );
            z = deltaz * ( i - 1 );
            
            % update material properties:
            % - if needed, k, S can be updated as function of r, z, Told
 
            % handle boundaries:
                if     i ==  1 % bottom
                C( i, j, ii ) = C0;
                        
                elseif i == Nz % top
                C( i, j , ii ) =  Cold( i-1, j ); 
                        
                elseif j ==  1 % left 
                C( i, j, ii ) =  Cold( i, j+1 );
                            
                        
                elseif j == Nr % right
                C( i, j, ii ) = ( Cold( i, j-1, ii ) + ( kw  * deltar / D ) * Cw  ) / ...
                            ( 1              + ( kw  * deltar / D ) );

                else % we are inside the domain:

                % coefficients:               
                Cleft   = -(D/(2*r*deltar)) - (D/(deltar^2));
                Cright  =  (D/(2*R*deltar)) - (D/(deltar^2));
                Cbottom = -(ax/(2*deltaz))  - (D/(deltaz^2));
                Ctop    =  (ax/(2*deltaz))  - (D/(deltaz^2));
                Ccenter =  Cbottom + Ctop +  Cleft  + Cright;
                
                % update solution:
                C( i, j, ii ) = ( Cleft   * Cold( i  , j-1, ii ) + ...
                              Cright  * Cold( i  , j+1, ii ) + ...
                              Cbottom * Cold( i-1, j, ii   ) + ...
                              Ctop    * Cold( i+1, j, ii   ) + ...
                              0 ) / Ccenter;
                          
                end
            end
        
        end
    
        % update error:
        error = norm( C(:,:,ii) - Cold(:,:,ii) );
    
        % display convergence:
        %fprintf('nite = %i, error = %2.2e \n ', error );

        % store the solution from this iteration:
        Cold(:,:,ii) = C(:,:,ii);

   end
end

%% Calculating Cavg & CO2 Absorbed %%

%extracting C values at top of pipe
C_TOP = C(N,:,:);

%Cavg calculation
syms r
C_avg = zeros(50,1);

for i = 1:length(C_TOP)

    C_avg(i) = ((1/(pi*R^2))*int(2*pi*r,r,0,R))*sum(C_TOP(:,:,i))/N;

end

C_avg

Answer 1

[Raag]

Hi Kyle,

As per my understanding the issue arises from incorrect indexing of the solution arrays with respect to the stack length variable in the loop.

• The result is the code repeatedly overwrites the same memory or uses the same values for each case, so all outputs are identical.
• To ensure that each stack length is treated independently and results are stored correctly, you should always index your solution arrays with the loop variable for stack length.
• When updating or referencing these arrays, use the third dimension for stack length.
• When averaging at the outlet, use the correct slice for each length.

This is how it can be done by replacing your main loop and averaging section with the following:

for ii = 1:num_L 
    % ... (your grid setup and boundary conditions) ... 
    C(:,:,ii) = 0;           % Solution for this stack length 
    Cold = zeros(Nz,Nr);     % Temporary storage 
 
    error = 1e8; 
    while error > tolerance 
        % ... (your iterative update, always use C(:,:,ii) and Cold) ... 
        % For example: 
        % C(i,j,ii) = ... (update using Cold, always with ii) 
        error = norm(C(:,:,ii) - Cold, 'fro'); 
        Cold = C(:,:,ii); 
    end 
end 
 
% Area-weighted average at the outlet for each length 
r_vec = linspace(0, R, Nr); dr = r_vec(2)-r_vec(1); 
C_avg = zeros(num_L,1); 
for ii = 1:num_L 
    C_top = C(Nz,:,ii); 
    C_avg(ii) = sum(2*pi*r_vec .* C_top)*dr / (pi*R^2); 
end 

This ensures each stack length is computed and stored independently, and your results will now vary as expected as shown in the output below:

For a better understanding of the above solution, refer to the following MATLAB documentations:

1. Array Indexing: https://www.mathworks.com/help/releases/R2021a/matlab/math/array-indexing.html

2. Partial differential equations: https://www.mathworks.com/help/releases/R2021a/matlab/math/partial-differential-equations.html