Let wave function $\Psi$ be defined on domain $D \in \mathbb{R}^n$. The Neumann condition $\frac{\partial \Psi} {\partial {\bf n}} = 0$ on the boundary $\partial D$ has a simple interpretation in terms of the probability current of $\Psi$. For $\Delta \Psi = i \partial\Psi/\partial t$ (although it's usually taken as $i \partial\Psi/\partial t = - \Delta \Psi$), the probability current at an arbitrary point ${\bf x} \in \mathbb{R}^n$ is
$$
{\bf j}({\bf x}) = i [ \Psi^*({\bf x}){\bf \nabla}\Psi({\bf x}) - \Psi({\bf x}){\bf \nabla}\Psi^*({\bf x}) ]
$$
and the normal current on $\partial D$ reads
$$
{\bf n} \cdot {\bf j} = i\; [ \Psi^* \frac{\partial \Psi}{\partial {\bf n}} - \Psi \frac{\partial \Psi^*}{\partial {\bf n}} ]
$$
(has the wrong sign, I know, but I accounted for OP's form of the Sch.eq. as $\Delta \Psi = i \partial\Psi/\partial t$).
Setting $\frac{\partial \Psi} {\partial {\bf n}} = 0$ amounts to ${\bf n} \cdot {\bf j} = 0$ everywhere on $\partial D$, thus confining the corresponding system within $D$ without an infinite potential well, as under Dirichlet conditions ($\Psi = 0$ on $\partial D$). This is the case of perfect reflection on $\partial D$.
There is a mention of this in Section 5.2 of Visual Quantum Mechanics:
Selected Topics with Computer-Generated Animations of Quantum-Mechanical Phenomena, by Bernd Thaller (Springer, 2000); Google Books link.
As for applications, one answer to another post, Can we impose a boundary condition on the derivative of the wavefunction through the physical assumptions?, pointed to the use of Neumann conditions in R-matrix scattering theory.
$$
$$
Clarification following @arivero's observation on conditions necessary to trap the system within domain D:
We can say that the system described by $\Psi$ is trapped in domain D if the total probability to locate it in D, $P_D$, is conserved in time: $dP_D/dt = 0$. In this case, if the system is initially located within D, such that $\Psi({\bf x},t=0) = 0$ for all ${\bf x} \notin D$ and $P_D(t=0) = 1$, then it will remain in D at all $t > 0$, since $P_D(t) = P_D(t=0) = 1$. If initially $P_D(t=0) < 1$ (the system has a nonzero probability to be located outside D), then we still have $P_D(t) = P_D(t=0) < 1$.
Conservation of $P_D$ is equivalent to a condition of null total probability current through the boundary $\partial D$. Note that it is not necessary to require null probability current at every point of $\partial D$, but only null total probability current through $\partial D$.
The difference can be understood in terms of path amplitudes (path-integral representation). In the former case, the amplitude $\Psi({\bf x}_1 t_1, {\bf x}_2 t_2; {\bf x} t)$ that the system "goes" from point ${\bf x}_1 \in D$ at time $t_1$ to point ${\bf x}_2 \notin D$ at time $t_2 > t_1$ while passing through point ${\bf x} \in \partial D$ at some time $t$, $t_1 < t < t_2$, is nonzero $\forall {\bf x} \in \partial D$: $\Psi({\bf x_1} t_1, {\bf x_2} t_2; {\bf x} t)≠0$. If however we demand null probability current at every point of $\partial D$, then $\Psi({\bf x}_1 t_1, {\bf x}_2 t_2; {\bf x} t)=0$, $\forall {\bf x} \in \partial D$.
In other words, null total probability current on $\partial D$ enforces weak trapping in the sense that overall $P_D(t) = $ const. and "cross-over events" across the boundary balance out. Null local probability current at every point of $\partial D$, ${\bf n} \cdot {\bf j} = 0$, corresponds to strong trapping in the sense that the system "does not cross" at any point ${\bf x} \in \partial D$. Imposing the strong trapping condition is equivalent to requiring that the weak trapping condition be satisfied by any wave function $\Psi$, as opposed to one selected $\Psi$. In this case the system is essentially confined within D at all times. Incidentally, the strong trapping condition follows from the requirement that the restriction of the system Hamiltonian on domain D remain self-adjoint.
Derivation of probability current conditions:
The free Schroedinger equation for $\Psi$, $i\partial\Psi/\partial t = \Delta\Psi$ as above (OP's choice of sign), implies local conservation of the probability density $\rho({\bf x,t}) = \Psi({\bf x},t)\Psi^*({\bf x},t)$:
$$
\frac{\partial \rho({\bf x},t)}{\partial t} + {\bf \nabla}\cdot {\bf j}({\bf x},t) = 0
$$
Integrating this over domain D yields
$$
\int_D dV\;\frac{\partial {\rho({\bf x},t)}}{\partial t} + \oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = \frac{d}{dt}\int_DdV\;{\rho({\bf x},t)} + \oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = 0
$$
which after denoting $P_D = \int_DdV\;{\rho({\bf x},t)}$ becomes
$$
\frac{dP_D}{dt} + \oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = 0
$$
Imposing $dP_D/dt = 0$ necessarily means $\oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = 0$. Note that $\oint_{\partial D} dS\;{{\bf n}\cdot{\bf j}} = 0$ does not require ${{\bf n}\cdot{\bf j}} = 0$ at every point on $\partial D$, whereas ${{\bf n}\cdot{\bf j}} = 0$ does imply $\oint_{\partial D} dS\;{{\bf n}\cdot{\bf j}} = 0$ and $dP_D/dt = 0$.
Self-adjoint restriction of the free Hamiltonian on domain D:
A self-adjoint restriction of $H\Psi = \Delta \Psi$ on D requires that $
\int_D dV\;\Phi^* (\Delta\Psi) = \int_D dV\;(\Delta\Phi^*) \Psi$ or $
\int_D dV\;[\Phi^* (\Delta\Psi) - (\Delta\Phi)^* \Psi] = 0$.
Use $\Phi^* (\Delta\Psi) = {\bf \nabla}\cdot(\Phi^* {\bf \nabla}\Psi) - {\bf \nabla}\Phi^* \cdot {\bf \nabla}\Psi$ and Green's theorem to obtain
$$
\int_D dV\;[\Phi^* (\Delta\Psi) - (\Delta\Phi)^* \Psi] = \oint_{\partial D} dS\;[\Phi^*\frac{d\Psi}{d{\bf n}} - \Psi\frac{d\Phi^*}{d{\bf n}}] = 0
$$
If the last condition above is to be satisfied by arbitrary $\Phi$, $\Psi$, it must hold locally:
$$
\Phi^* \frac{d \Psi}{d{\bf n}} - \Psi \frac{d \Phi^*}{d{\bf n}} = 0
$$
This means $\frac{1}{\Psi}\frac{d\Psi}{d{\bf n}} = \frac{1}{\Phi^*}\frac{d\Phi^*}{d{\bf n}} = a({\bf x}), \forall {\bf x} \in \partial D$. The case $\Phi = \Psi$ shows that $a({\bf x}) = a^*({\bf x})$. Therefore the restriction of $H$ on $D$ is self-adjoint if and only if wave functions in its support satisfy a boundary condition
$$
\frac{d\Psi}{d{\bf n}} = a({\bf x})\Psi
$$
for some given real-valued function $a({\bf x})$. In particular, this means every wave function $\Psi$ also satisfies the strong trapping condition ${\bf n} \cdot {\bf j} = 0$.
Finally note that the strong trapping condition means $\Psi^*\frac{d\Psi}{d{\bf n}} - \Psi\frac{d\Psi^*}{d{\bf n}} = 0$, but does not imply that $\Phi^*\frac{d\Psi}{d{\bf n}} - \Psi\frac{d\Phi^*}{d{\bf n}} = 0$ for arbitrary $\Psi$, $\Phi$. If we consider the latter expression as the matrix element of a "local current operator" ${\bf n} \cdot {\bf \hat j}$, then the strong trapping condition requires that diagonal elements of ${\bf n} \cdot {\bf \hat j}$ are 0, whereas self-adjointness requires that the entire space of wave functions is in the kernel of each ${\bf n} \cdot {\bf \hat j}$.
Best Answer
Here is a simple variant on the proof in the OP that occurred to me after posting the bounty. I am once again considering two solutions $\phi'$ and $\phi''$, with $\vec{D} = \epsilon \vec{E} = -\epsilon \nabla \phi$.
The divergence theorem and a little rearranging nets us $\int_{\Gamma} dV \nabla(\phi' - \phi'') (\vec{D}'-\vec{D}'') = \int_{\partial \Gamma}(\phi' - \phi'') (\vec{D}'-\vec{D}'')\cdot \vec{dA} - \int_{\Gamma} dV (\phi' - \phi'') \nabla \cdot (\vec{D}'-\vec{D}'')$
Then, once again, for Dirichlet boundary conditions, $\phi' - \phi''$ is zero on the boundary, while for von Neumann boundary conditions, $\vec{D}' - \vec{D}''$ is zero on the boundary (since it's proportional to $0$). Thus, for these boundary conditions or an appropriate mixture, the boundary term vanishes.
Furthermore, $\nabla \cdot \vec{D}' =\nabla \cdot \vec{D}'' = \rho_{\rm free,\, specified}$, so the second term on the right hand side also vanishes. Thus, $$\int_{\Gamma} dV \nabla(\phi' - \phi'') (\vec{D}'-\vec{D}'') = 0$$ which means $$\int_{\Gamma} dV \epsilon |\nabla(\phi' - \phi'')|^2 = 0$$
If $\epsilon$ does not change sign, then $\nabla(\phi' - \phi'')=0$, proving uniqueness of the electric field.