Introduction
This code is useful in the field of Environmental and Hydrogeology Sciences. Where most groundwater problems appear in most cases such as calculating the groundwater budget under the stresses of discharge/recharge (pumped or injected wells).
Background
Finite difference approximation is a fundamental numerical math method that helps to solve the problem of complex equations that function natural physical problem. Especially in the problem of transient condition (time dependent function), and there are three general widely used finite difference schemes.
1. Explicit Differentiation
h(i,j,k+1) = rh(i-1,j,k)+(1-4r)h(i,j,k)+rh(i+1,j,k)+rh(i,j-1,k)+rh(i,j+1,k)
2. Implicit
h(i,j,k) = rh(i-1,j,k+1)+(1+4r)h(i,j,k+1)+rh(i+1,j,k+1)+rh(i,j-1,k+1)+rh(i,j+1,k+1)
3. Crank Nicholson
rh(i-1,j,k)+(2-4r)h(i,j,k)+rh(i+1,j,k)+rh(i,j-1,k)+rh(i,j+1,k) =
rh(i-1,j,k+1)+(2+4r)h(i,j,k+1)+rh(i+1,j,k+1)+rh(i,j-1,k+1)+rh(i,j+1,k+1)
And finally, combining all three methods by adding Lamda coefficient (code below):
h(i,j,k+1) = [h(i,j,k)+r*lamda*[h(i-1,j,k)-4h(i,j,k)+h(i+1,j,k)+
h(i,j-1,k)+h(i,j+1,k)]+r(1-lamda)*[h(i-1,j,k+1)+h(i+1,j,k+1)+h(i,j-1,k+1)+h(i,j+1,k+1)]]/(1+4r(1-lamda)
where:
- r = D*DT/DX^2
- D = T/S
- T: is the transmissivity of the aquifer (m2/day)
- S: aquifer sorage coefficient
- DX: change in distance in X-direction
- DT: time increment
- lamda:
- if = 0 ----> Implicit, if = 0.5 ----> Crank Nicholson, if = 1 ----> Explicit
Using the Code
In this code, we will solve the groundwater head of confined aquifer distribution in which a sudden drop occurs to groundwater at initial condition was 8.8 m distributed in all the aquifer domain. The water flow from upstream boundary with constant head of 8.8 m towards the downstream boundary 3.2 m.
The figure below illustrates the conceptual model of the aquifer:
In the code below, Gauss Seidel Iteration method was applied to solve the finite difference scheme.
double T = 2; double S = 0.0003; double DT = 0.5; double DX = 100; double lamda = 0; double D = T / S;
double H0 = 8.8; double H1 = 8.8; double H2 = 3.2; int IT = 100; double ERR = 1; double TOL = 0.001; double r = (D * DT) / System.Math.Pow(DX, 2);
double[, ,] H = new double[10, 10, 18];
for (int i = 0; i <= 9; i++)
{
for (int j = 0; j <= 9; j++)
{
int k = 0;
H[i, j, k] = H0;
}
}
for (int i = 0; i <= 9; i++)
{
int j = 0;
for (int k = 1; k <= 17; k++)
{
H[i, j, k] = H1;
}
}
for (int i = 0; i <= 9; i++)
{
int j = 9;
for (int k = 1; k <= 17; k++)
{
H[i, j, k] = H2;
}
}
while (ERR > TOL)
for (int n = 0; n < IT; ++n)
{
for (int i = 1; i <= 8; i++)
{
for (int j = 1; j <= 8; j++)
{
for (int k = 0; k <= 16; k++)
{
double OLD = H[i, j, k + 1];
H[i, j, k + 1] = (H[i, j, k] + (r * lamda) * (H[i - 1, j, k] -
(4 * H[i, j, k]) +
H[i + 1, j, k] + H[i, j - 1, k] + H[i, j + 1, k]) +
(r * (1 - lamda)) * (H[i - 1, j, k + 1] + H[i + 1, j, k + 1] +
H[i, j - 1, k + 1] +
H[i, j + 1, k + 1])) / (1 + (4 * r * (1 - lamda)));
H[0, j, k + 1] = H[2, j, k + 1];
H[9, j, k + 1] = H[7, j, k + 1];
ERR = System.Math.Abs(H[i, j, k + 1] - OLD);
}
}
}
}
Console.WriteLine("Groundwater Head values in (m)");
Console.WriteLine(IT.ToString());
Console.WriteLine(ERR.ToString());
for (int i = 0; i <= 9; i++)
{
for (int j = 0; j <= 9; j++)
{
int k = 17;
Console.Write(string.Format("{0} ", H[i, j, k]));
}
Console.Write(Environment.NewLine + Environment.NewLine);
}
Console.ReadLine();
}
}
}
Output
Head distribution at the end of time (t = 18
), with number of iterations = 100
tolerance = 0.002
Points of Interest
This code was implemented with C# language using console application to solve the PDE in the environmental sciences with the application of FDM in two-dimensions.
References
- Mary P Anderson. Introduction to Groundwater Modeling: Finite Difference Finite Element
- Amr El-Feki. Prof in Environmental Sciences K A U
History
- 4th October, 2018: Initial version
This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.