## 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)

h(:,end)

toc(t)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) 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) finite difference method, lapplace equation, central difference method, numerical analysis MATLAB Answers — New Questions