Mineral dissolution and solute transport around a solid

Question

I am trying to simulate solute transport of acid (HCl) and consequent mineral dissolution around a grain (calcite).

The governing equation for transport is the advection-diffusion equation, given as:

where C is the concentration, v is the flow velocity, D is the diffusion coefficient. At a constant flow velocity, the above is discretized using an FTCS formulation as:

The dissolution occur at the fluid/solid interface – and can be modelled as a Robin type Bc, which is described as:

where r is the mineral rate of dissolution, ε is the stoichiometric coefficient of species in the dissolution process (this is equal to 1 in the reaction considered), k_(H+) is the rate constant, γ_(H+) is the activity coefficient, and c_(H+) is the molar concentration of reacting specie (in this case Hydrogen or chloride ion).

I have discretized Eq. (4) as:

This is my code so far:

clc; clear; close all;


%%%%%%%%%% Specify inputs
CircleDiam = 40;        %in pixels
Co    = 0;          %initial conc. in the domain [mol/L]
Cin   = 0.01;       %conc. of injected fluid [mol/L]

Lx    = 0.1/100;        %Length of domain [m]
nx    = 200;        %spatial gridpoints in x
dx    = Lx/(nx-1);      %Length step size [m]
Ly    = 0.05/100;       %Width of domain [m]
ny    = 100;        %spatial gridpoints in y
dy    = Ly/(ny-1);      %Length step size [m]

T     = 5/(3600*24);    %Simulation time [days]
nt    = 8000;           %shifts
dt    = T/nt;       %Time step [days]

% Flow
u = 103.68;         %Velocity in x direction [m/day]
v = 0;              %Velocity in y direction [m/day]
De = 8.64e-04;      %Dispersion coeff. [m2/day]
betaX  = u*dt/dx;
betaY  = v*dt/dy;
gammaX = De*dt/(dx^2);
gammaY = De*dt/(dy^2);

%%%%%%%%%% Create image with solid object
radius  = CircleDiam/2;

% obtain full output grids from grid vectors
[Colgrids, Rowgrids] = meshgrid(1:nx, 1:ny);

% create a logical mask for the circle by specifying the center and diameter of the circle.
centerX = 0.5 + (nx/2);
centerY = 0.5 + (ny/2);

% obtain image from: ( (y-y0)^2 + (x-x0)^2 )  <= r^2, where (y0,x0) is the centre point of circle
SolidImg = (Rowgrids - centerY).^2 + (Colgrids - centerX).^2 <= radius.^2;

% change from logical to numeric labels. Also, transpose matrix to conform with the conc. matrix
P   = double(SolidImg');
%figure, imshow(~SolidImg, [], 'InitialMagnification','fit'); box on;


%%%%%%%%%% Find the solid boundaries
BW = im2bw(P);
BW_filled  = imfill(BW,'holes');
boundaries = bwboundaries(BW_filled);

for k = 1:size(boundaries, 1)
    bdr_matrix = boundaries{k};
end

row = bdr_matrix(:,1);
col = bdr_matrix(:,2);

% Gridblocks
x = 0:dx:Lx;
y = 0:dy:Ly;
t = 0:dt:T;
[X,Y] = meshgrid(x, y);

% specify initial conditions
C        = zeros(nx, ny, nt+1);
C(:,:,1) =  Co;              %Initial condition

% specify Reaction parameters
kh     = 77004.081;     % [mol/m2/day]
gammaH = 1;             % [L/mol]


%iterate finite difference equations
for k = 1:nt
        for j = 2:ny-1
    for i = 2:nx-1
        if P(j,i)==1    %Solid pixels
            C(j,i,k+1) = C(j,i,k);

            % Apply the robin bc at the fluid/solid interface
            if ismember(P(j,i), P(row, col))
                C(j,i,k+1) = C(j,i,k)...
                    - 1000*0.5*De*(C(j,i+1,k)-C(j,i-1,k)/dx)/(kh*gammaH)...
                    - 1000*0.5*De*(C(j+1,i,k)-C(j-1,i,k)/dy)/(kh*gammaH); % multiply by 1000L to make dimensionally consistent with mol/L
            end
        else
            C(j,i,k+1) = C(j,i-1,k)*(betaX/2+gammaX) + C(j,i+1,k)*(gammaX-betaX/2)...
                + C(j-1,i,k)*(betaY/2+gammaY) + C(j+1,i,k)*(gammaY-betaY/2)...
                + C(j,i,k)*(1-2*gammaX-2*gammaY);
        end
    end
end

% Insert boundary conditions   
C(1,:,k+1)   = C(2,:,k+1);      % south
C(end,:,k+1) = C(end-1,:,k+1);  % north
C(:,1,k+1)   = Cin;             % west
C(:,end,k+1) = C(:,end-1,k+1);  % east
end

My result (see image) does not look like what I would expect (see second image). Ideally, I expect to see something like the above image, and so I am wondering where I might be getting it wrong. I have assumed that both velocity and diffusion coefficient are zero at the solid cells, and I have also prescribed a robin type bc at the fluid/solid interface, however, I get the same outcome when I run the code without eq. (5). I want to believe that I am not doing this the right way. I am not completely sure if I have discretized Eq. (5) the right way, and I am wary of how I have implemented it (i.e., eq. 5) in the code. I thus would very much appreciate if someone can help me out here or point me in the right direction.

PS: The axes are the same in both image. The difference is that mine is in pixels while the other is in cm.

Thank you in anticipation.

EDIT

Thank you so much for your help @Sthavishtha Bhopalam . I am indeed very grateful.

I do have some questions/comments/clarifications:

1. I obtained the velocity field from LBM. It’s a pity I didn’t mention this initially. I have now included the appropriate lines of code (without the LBM though) as shown below - and I'd appreciate if you can confirm the correctness.

The following comes just before the loop:

 % calculate flow velocity field
    [u,v]  = FlowVelocitySolver(nx,ny,1000);  % iter = 1000;

And then the values of u and v at every point in the grid is utilized in the advection term.

2. I have run the code with your amendments – and I got the expected shape as in Molins et al. (see the new figure), however:

• I see you changed the inlet velocity from 103.68 m/day (which is 0.12 cm/s) to 10.368 m/day (which is 0.012 cm/s), the former being the same with Molins et al.

• I noticed that numerical instability occurs at the fluid-solid interface when I use 103.68 m/day. This happens even after I adjusted the time step value to match the stability criteria. Is this the reason you changed to 10.368 m/day?

3. I was looking to use other kinetic formulations. For example, see pages 68-69 of: https://biblio.ugent.be/publication/8071415/file/8071428.pdf

Would it be right to discretize Eqs. 4.10, given by:

$ -D\textbf n.\nabla c = \frac{V}{A_{s}}k_{het}(c – c_{eq}) $

as:

Many thanks once again.

Answer 1

Firstly, there are some mistakes in your discretized formulation and the approach you employ to solve this problem:

• The flow velocity in the advection diffusion Eq. (1) is to be computed from the fluid continuity (COM) and momentum equations (COLM). In other words, you cannot assume the flow velocity to be constant in the entire domain, i.e., you employ constant variables u and v in the code. The correct result (which you extract from ref. [1]) is intuitively identical to the flow over a two-dimensional cylinder at a low Reynolds number. So, this necessitates the need to solve COM and COLM.
• To ensure numerical stability, the advective part of Eq. (1) has to be discretized by an upwind scheme. So, the discrete form would be (for the purpose of illustration here, I assume $u > 0, v > 0$)

$$ \frac{c_{i,j}^{n+1} - c_{i,j}^{n}}{\Delta t} = D \left(\frac{c_{i+1,j}^n - 2c_{i,j}^n + c_{i-1,j}^n}{(\Delta x)^2} + \frac{c_{i,j+1}^n - 2c_{i,j}^n + c_{i,j-1}^n}{(\Delta y)^2}\right) - \frac{(uC)_{i,j}^n - (uC)_{i-1,j}^n}{\Delta x} - \frac{(vC)_{i,j}^n - (vC)_{i,j-1}^n}{\Delta y} $$

• After obtaining the discrete equation, you must appropriately choose the input parameters to ensure the satisfaction of the CFL criterion, i.e., positivity of the coefficient of $c_{i,j}^n$.

• The discrete form Eq. (5) of the Robin BC Eq. (4) needs to be changed. The LHS of the discrete form would be (I just show forward difference here, but you could also use central differencing)

$$ -D\left(n_x \frac{C_{i,j}^n - C_{i-1,j}^n}{\Delta x} + n_y \frac{C_{i,j}^n - C_{i,j-1}^n}{\Delta y}\right) = RHS $$

where $n_x, n_y$ are the unit normal vectors to the fluid-solid interface. You may have to double check the direction of these vectors from appropriate references. Also, I think that the term $c_{H+}$ in the RHS is the concentration at the fluid-solid interface, which would not be equivalent to $c_{i,j}^{n+1} - c_{i,j}^{n}$.

• I have made some of the aforementioned changes to your script, except of solving for COM and COLM. As you may see, the solution still partially resembles your plots, primarily because of neglecting these equations.

1. Molins, Sergi; Soulaine, Cyprien; Prasianakis, Nikolaos I.; Abbasi, Aida; Poncet, Philippe; Ladd, Anthony J. C.; Starchenko, Vitalii; Roman, Sophie; Trebotich, David; Tchelepi, Hamdi A.; Steefel, Carl I., Simulation of mineral dissolution at the pore scale with evolving fluid-solid interfaces: review of approaches and benchmark problem set, ZBL07420582.

    clc; clear; close all;
    % parameters slightly modified as per the paper
    %%%%%%%%%% Specify inputs
    CircleDiam = 40;    %in pixels
    Co    = 0;          %initial conc. in the domain [mol/L]
    Cin   = 0.01;       %conc. of injected fluid [mol/L]
    
    Lx    = 0.1/100;    %Length of domain [m]
    nx    = 200;        %spatial gridpoints in x
    dx    = Lx/(nx-1);  %Length step size [m]
    Ly    = 0.05/100;   %Width of domain [m]
    ny    = 100;        %spatial gridpoints in y
    dy    = Ly/(ny-1);  %Length step size [m]
    
    T     = 500/(3600*24);%Simulation time [days]
    nt    = 5E4;        %shifts
    dt    = T/nt;       %Time step [days]
    
    % Flow
    % u = 103.68;       %Velocity in x direction [m/day]
    u = 10.368;
    v = 0;              %Velocity in y direction [m/day]
    De = 8.64e-06;      %Dispersion coeff. [m2/day]
    betaX  = u*dt/dx;
    betaY  = v*dt/dy;
    gammaX = De*dt/(dx^2);
    gammaY = De*dt/(dy^2);
    
    %%%%%%%%%% Create image with solid object
    radius  = CircleDiam/2;
    
    % obtain full output grids from grid vectors
    [Colgrids, Rowgrids] = meshgrid(1:nx, 1:ny);
    
    % create a logical mask for the circle by specifying the center and diameter of the circle.
    centerX = 0.3 + (nx/2);
    centerY = 0.3 + (ny/2);
    
    % obtain image from: ( (y-y0)^2 + (x-x0)^2 )  <= r^2, where (y0,x0) is the centre point of circle
    SolidImg = (Rowgrids - centerY).^2 + (Colgrids - centerX).^2 <= radius.^2;
    
    % change from logical to numeric labels. Also, transpose matrix to conform with the conc. matrix
    P   = double(SolidImg');
    %figure, imshow(~SolidImg, [], 'InitialMagnification','fit'); box on;
    
    
    %%%%%%%%%% Find the solid boundaries
    BW = im2bw(P);
    BW_filled  = imfill(BW,'holes');
    boundaries = bwboundaries(BW_filled);
    
    for k = 1:size(boundaries, 1)
        bdr_matrix = boundaries{k};
    end
    
    % gxn and gyn have been accidentally swapped
    [gmag,gdir] = imgradient(BW_filled);
    [gx,gy] = imgradientxy(BW_filled);
    gxn = gx;
    gyn = gy;
    
    % normalized gradient
    for i = 1:nx
        for j = 1:ny
            gxn(i,j) = gx(i,j)/gmag(i,j);
            gyn(i,j) = gy(i,j)/gmag(i,j);
        end
    end
    
    row = bdr_matrix(:,1);
    col = bdr_matrix(:,2);
    
    % Gridblocks
    x = 0:dx:Lx;
    y = 0:dy:Ly;
    t = 0:dt:T;
    [X,Y] = meshgrid(x, y);
    
    % specify initial conditions
    C        = zeros(nx, ny);
    C(:,:)   = Co;              %Initial condition
    Cnew     = C;
    
    % specify Reaction parameters
    kh     = 77004.081;     % [mol/m2/day]
    gammaH = 1;             % [L/mol]
    
    %iterate finite difference equations
    for k = 1:nt
            for j = 2:ny-1
            for i = 2:nx-1
                if P(i,j)==1    %Solid pixels
                    Cnew(i,j) = 0.;
    
                    % Apply the robin bc at the fluid/solid interface
                    for l = 1:size(row,1)
                            if (i == row(l)) && (j == col(l))                           
                            temp = ...
                                - 1000*De*.5*gyn(i,j)*(C(i+1,j)-C(i-1,j)/dx)/(kh*gammaH)...
                                - 1000*De*.5*gxn(i,j)*(C(i,j+1)-C(i,j-1)/dy)/(kh*gammaH); % multiply by 1000L to make dimensionally consistent with mol/L
                            Cnew(i,j) = temp;
                            end
                    end
                else
                    Cnew(i,j) = C(i-1,j)*(betaX/2+gammaX) + C(i+1,j)*(gammaX-betaX/2)...
                            + C(i,j-1)*(betaY/2+gammaY) + C(i,j+1)*(gammaY-betaY/2)...
                            + C(i,j)*(1-2*gammaX-2*gammaY);
                end
            end
            end
    
        % Insert boundary conditions   
        Cnew(:,1)   = Cnew(:,2);      % south
        Cnew(:,end) = Cnew(:,end-1);  % north
        Cnew(1,:)   = Cin;            % west
        Cnew(end,:) = Cnew(end-1,:);  % east
    
        k = k + 1;    
        C = Cnew;
    end