Finite Difference Method - Central Difference Method for Groundwater Flow Problem

Hi all!

I am a student who is trying to make a MATLAB code to solve a simple steady-state groundwater flow problem using the finite difference method-central difference method. The boundary is a square boundary with a dimension of 10 cm x 10 cm. The uppermost and lowermost boundaries are impermeable. The total constant head at the leftmost and rightmost boundaries are 10 cm and 5 cm, respectively. However, the results are somehow not correct because the total head at the leftmost and rightmost boundaries are not equal to 10 cm and 5 cm, respectively. I suspect that there is something wrong with my boundary condition and I did not know how to solve it. Your help is appreciated. Here is my code:

%This program is designated to calculate 2D seepage flow in square domain
%The case is for isotropic homogeneous condition.
clc;  clear all; close all;
t=tic;
%% Define the domain and soil properties
Lx = 10; %Length of the domain in the x-direction in cm
Ly = Lx; %Length of the domain in the y-direction in cm; the domain is square
kxx = 1; kyy =1; kxy = 0; kyx = kxy; %Soil properties for isotropic condition with kxy = kyx = 0 (in cm/s)
Nx = 11; %Number of grid points in x-direction
Ny = Nx; % Number of grid points in y-direction, the same as Nx

%Discretization
dx = Lx/(Nx-1); % Grid spacing in x-direction
dy = Ly/(Ny-1); % Grid spacing in y-direction

%% Head boundary conditions
h_left = 10; % Left boundary head (in cm)
h_right = 5; % Right boundary head (in cm)

%% Defining the the governing equation in the entire domain
Ix = speye(Nx); Iy = speye(Ny);
Kx = spdiags(ones(Nx,1)*[kxx -2*kxx kxx],-1:1,Nx,Nx);
Ky = spdiags(ones(Ny,1)*[kyy -2*kyy kyy],-1:1,Ny,Ny);
A = kron(Ky/dy^2,Ix)+kron(Iy,Kx/dx^2);
A = full(A);

%% Defining boundary conditions
% Initial Boundary Condition
bc = zeros(length(A),1);

%Applying the top impermeable boundary
A(1:Nx,1+Ny:2*Ny)=2*A(1:Nx,1+Ny:2*Ny);
%Applying the bottom impermeable boundary
A(Nx+(Nx-1)*(Ny-1):1:(Nx*Ny),(Nx-1)*(Ny-1):1:(Ny-1)+(Nx-1)*(Ny-1))=...
    2*A(Nx+(Nx-1)*(Ny-1):1:(Nx*Ny),(Nx-1)*(Ny-1):1:(Ny-1)+(Nx-1)*(Ny-1));
%Left boundary
A(1:Nx:Nx+(Nx-1)*(Ny-1),1)=2*A(1:Nx:Nx+(Nx-1)*(Ny-1),1);
%Right boundary
A(Nx:Nx:Nx+(Nx-1)*(Ny-1),Nx)=2*A(Nx:Nx:Nx+(Nx-1)*(Ny-1),Nx);

% Applying left boundary (Dirichlet BC)
A(1:Ny:Nx*Ny,1:Ny:Nx*Ny)=0;
bc(1:Ny:Ny*Ny,1)= h_left;
% Ayplying right boundary (Dirichlet BC)
A(Ny:Ny:Nx*Ny,Ny:Ny:Nx*Ny)=0;
bc(Nx:Nx:Nx*Ny,1)= h_right;

%% Solve the PDE using Direct Matrix Inversion
h = inv(A)*bc;
h = reshape(h, Nx, Ny)'; % Reshape solution to 2D grid

%% Plotting the solution
% Plotting the sparse matrix
figure; spy(A);
xlabel('Column index'); ylabel('Row index');
title('Sparse Coefficient Matrix (A)');

[x, y] = meshgrid(0:dx:Lx, 0:dy:Ly);
colormap(jet);

contourf(x, y, h, 20, 'LineColor', 'none');
colorbar;
xlabel('X (m)');
ylabel('Y (m)');
title('Groundwater Flow Problem');

h(:,1)

ans = 11x1

10.6250 10.6250 10.6250 10.6250 10.6250 10.6250 10.6250 10.6250 10.6250 10.6250

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

h(:,end)

ans = 11x1

4.3750 4.3750 4.3750 4.3750 4.3750 4.3750 4.3750 4.3750 4.3750 4.3750

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

toc(t)

Elapsed time is 2.200534 seconds.

Answer 1

[Torsten]

It would help if you included the mathematical description of your problem in order to compare with your coding (equation(s) to be solved, initial condition(s),boundary conditions).

A Dirichlet boundary condition is simply set as

h(:,1) = h_left
h(:,end) = h_right

Thus your coefficient matrix A must have a "1" for the coefficient related to the boundary point and h_left resp. h_right for the corresponding value of bc.

The solution of your problem is just a linear function f between x=0 and x=10 with f(0) = h_left and f(10) = h_right independent of the y-coordinate.

[Ignatius Tommy Pratama]

The following attachment shows the theoretical explanations. Basically, I want to replicate the problem in the attached PPT. But, somehow, I stuck at the boundary conditions. I really appreciate your help.

Thank you

[Torsten]

Do you know how the strange solution from slice 11 is produced ? Is it a time-dependent simulation ?

%% Define the domain and soil properties
Lx = 10; %Length of the domain in the x-direction in cm
Ly = Lx; %Length of the domain in the y-direction in cm; the domain is square
kxx = 2; kyy = 1.5; kxy = 1; kyx = kxy; %Soil properties for isotropic condition with kxy = kyx = 0 (in cm/s)
nx = 11; %Number of grid points in x-direction
ny = nx; % Number of grid points in y-direction, the same as Nx
%Discretization
dx = Lx/(nx-1); % Grid spacing in x-direction
dy = Ly/(ny-1); % Grid spacing in y-direction
x = 0:dx:Lx;
y = 0:dy:Ly;
%% Head boundary conditions
h_left = 10; % Left boundary head (in cm)
h_right = 5; % Right boundary head (in cm)

a = kxx/dx^2; 
c = kyy/dy^2;
b = -2*(a+c);
d = (kxy+kyx)/(4*dx*dy);
A = zeros(nx*ny,nx*ny);
B = zeros(nx*ny,1);
% Boundaries
% Boundary values at y = 0 (Impermeable boundary)
for ix = 2:nx-1
    k = ix;
    A(k,k) = b;
    A(k,k+1) = a;
    A(k,k-1) = a;
    A(k,k+nx) = 2*c;
    B(k) = 0.0;
end
% Boundary values at y = 1 (Impermeable boundary)
for ix = 2:nx-1
    k = nx*(ny-1) + ix;
    A(k,k) = b;
    A(k,k+1) = a;
    A(k,k-1) = a;
    A(k,k-nx) = 2*c;
    B(k) = 0.0;
end
% Boundary values at x = 0 (H = h_left)
for iy = 1:ny
    k = nx*(iy-1) + 1;
    A(k,k) = 1.0;
    B(k) = h_left;
end
% Boundary values at x = 1 (H = h_right)
for iy = 1:ny
    k = nx*iy;
    A(k,k) = 1.0;
    B(k) = h_right;
end

% Inner grid points
for iy = 2:ny-1
    for ix = 2:nx-1
        k = (iy-1)*nx + ix;
        A(k,k) = b;
        A(k,k+1) = a;
        A(k,k-1) = a;
        A(k,k+nx) = c;
        A(k,k-nx) = c;
        A(k,k+1+nx) = d;
        A(k,k-1+nx) = -d;
        A(k,k+1-nx) = -d;
        A(k,k-1-nx) = d;
    end
end
h = A\B;
%u
H = zeros(nx,ny);
for iy = 1:ny
    for ix = 1:nx
        k = (iy-1)*nx + ix;
        H(ix,iy) = h(k);
    end
end
H(1,:)

ans = 1x11

10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

H(end,:)

ans = 1x11

5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% Plotting the solution
% Plotting the sparse matrix
figure; spy(A);
xlabel('Column index'); ylabel('Row index');
title('Sparse Coefficient Matrix (A)');

[X,Y] = meshgrid(x,y);
colormap(jet);

contourf(X,Y,H.', 20, 'LineColor', 'none');
colorbar
xlabel('X (m)');
ylabel('Y (m)');
title('Groundwater Flow Problem');

[Ignatius Tommy Pratama]

Thank you very much for your help. Have a great day!