# Finite difference aproximation - Darcy law

I am solving following problem:

Filtration of water can be described in bi-dimensional case by $$- \partial_x(K(x,y) \partial_x u ) - \partial_y (K(x,y) \partial_y u ) = 0,$$ where $u$ - water level $K$ - hydraulic conductivity On uniform mesh of $N$, $K$ is defined as \begin{align} K(x,y) & = \begin{cases} (x -0.5)^2 + (y - 0.5)^2 - 0.04 & \geq 0\\ & <0 \end{cases}\\ \end{align}

\begin{align} - \partial_x K \partial_x u \ \ \text{is defined as} \ \ \frac{K(i+ 1/2) \partial_x u_{i+1/2} - K(i - 1/2) \partial_x u_{i - 1/2} }{h}, \end{align}

where $(i + 1/2)$ refers to mid-points of $(i + 1)th$ interval.

I have extended discretization into 2D case and obtained following system (not sure, if its correct) \begin{align} \begin{bmatrix} c_{i} & c_{i+1} & 0 & d_{i+1} & 0 & 0 & 0 & 0 \\ c_{i-1} & c_{i} & c_{i+1} & 0 & d_{i+1} & 0 & 0& 0 \\ 0 & c_{i-1} & c_{i} & c_{i+1} & 0 & d_{i+1}& 0 & 0\\ d_{i-1} & \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \ddots \\ 0 & 0 & 0 & \dots &\dots & & c_{i-1} & c_{i} \end{bmatrix} * \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ u_{N-1, 1} \\ \vdots \\ u_{N-1,N-1} \\ \end{bmatrix} = 0 \end{align} where coefficients are calculated as \begin{align} c_i &= + K(i -1/2,j) + K(i +1/2,j) + K(i ,j+1/2) + K(i,j-1/2)\\ c_{i + 1} &= - K(i+1/2,j) \\ c_{i- 1} &= - K(i - 1/2,j) \\ d_{i - 1} &= - K(i,j - 1/2) \\ d_{i + 1} &= - K(i,j +1/2) \\ \end{align}

System has to satisfy following BC: \begin{align} \partial_x u &= 0 \ \ &\text{on} \{0,1 \} \times (0,1)\\ u &= 10 \ \ &\text{on} \{0 \} \times (0,1) \\ u &= 0 \ \ &\text{on} \{1 \} \times (0,1) \end{align}

Can anybody give me some advices how to handle them ?

• What confuses you the most about implementing these boundary conditions?
– Paul
Oct 31 '15 at 16:17
• @Paul lot of staff... what I was thinking was: for Dirichlet = 0, I can just skip it from system, for neuman I am not sure what to do with them and for Dirichlet = 10 : after solving system assign them, so in this way I could preserve symmetry of system... Oct 31 '15 at 16:23
• Take a look at this answer on a similar question, specifically the section Adding a convective boundary condition. If you set $h=0$ you have a zero-gradient condition. Maybe it is helpful. Oct 31 '15 at 23:21