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
