3
$\begingroup$

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:

enter image description here

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:

enter image description here

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

enter image description here

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:

enter image description here

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

enter image description here

enter image description here

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.

  1. 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?

  1. 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: enter image description here

Many thanks once again.

enter image description here

$\endgroup$
2
  • $\begingroup$ I tried running this in Octave, an open source Matlab clone, but nothing came out. And I don't want to convert it to Python by hand to poke at it. I'd like to help, but I don't want to buy Matlab. Your choice of programming language has made it difficult for me to help you :-( $\endgroup$
    – Richard
    Dec 2, 2021 at 7:16
  • 1
    $\begingroup$ @Richard, my bad I am still a python learner at the moment. Thank you for even trying to help. Cheers! $\endgroup$
    – oma11
    Dec 2, 2021 at 11:05

1 Answer 1

2
$\begingroup$

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
$\endgroup$
6
  • $\begingroup$ you should post the questions/comments related to my answer here, instead of editing in the original post. Reg. the flow velocity, have you tried changing the diffusion coefficient to $8.64e-5$, as used in the paper. you need to do the same for other parameters you employ in your lbm code. $\endgroup$ Dec 14, 2021 at 23:54
  • $\begingroup$ also, I still spot some mistakes in the changes you made in the code. $(uC)_{i-1,j}$ is not equivalent to $u_{i,j}C_{i-1,j}$, as $u$ is varying (obtained from your flow solver). $\endgroup$ Dec 15, 2021 at 0:00
  • $\begingroup$ reg. eqn. 4.10: it looks like $c - c_{eq}$ should actually be $c_{\Gamma} - c_{eq}$, where $c_{\Gamma}$ denotes the concentration at the fluid-solid interface. and $c_{\Gamma}$ should be somehow linked to the solid concentration; but i am not very sure. $\endgroup$ Dec 15, 2021 at 0:29
  • $\begingroup$ Thanks again for your response. You're right, I shouldn't have added anything to your code. So, I have now removed my additions as appropriate. Furthermore, I didn't realize you changed De from Molins. In any case, I thought the original value should be 8.64e-4 ? Or I'm I missing something? $\endgroup$
    – oma11
    Dec 15, 2021 at 13:01
  • $\begingroup$ Thank you also for your comment regarding eqn. 4.10. As pertaining the expansion of $(uC)_{i-1,j}$, please what should be the correct thing? $\endgroup$
    – oma11
    Dec 15, 2021 at 13:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.