Finite Difference Method - Central Difference Method for Groundwater Flow Problem
Show older comments
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)
h(:,end)
toc(t)
Accepted Answer
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.
3 Comments
Ignatius Tommy Pratama
on 9 May 2024
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,:)
H(end,:)
%% 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
on 11 May 2024
More Answers (0)
Categories
Find more on Agriculture in Help Center and File Exchange
on 8 May 2024
on 11 May 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
