This will only be a partial answer to your question, due to lack of time.

As you effectively noted, any $f\in C^2(\overline{D})$ with compact support and such that $(\nu, \nabla f) = 0 \ \mathcal{H}^{d-1}$-a.e. on $\partial D$ will belong to $D(L)$. The problem with your argument is that, given a differentiability point $p \in \partial D$ and a general $U \subseteq \mathbb{R}^d$ and $\psi : B(1) \to U$ as you describe, it is not true that $$(\nu, \nabla f) = 0 \ \text{at} \ p \Leftrightarrow \partial_d(f \circ \psi)\big|_0 = 0 .$$ To see why, suppose you had a different pair $(U',\psi')$ satisfying the same requirements. Then, by the chain rule, a.e. on $\psi'^{-1}(U \cap U')$ it holds that $$\partial_d(f \circ \psi')\big|_{x} = \partial_d(f \circ \psi \circ \psi^{-1} \circ \psi')\big|_{x} = \sum_{i=1}^d \partial_i(f \circ \psi)\big|_{y(x)}\partial_d y_i\big|_{x}$$ where $y := \psi^{-1} \circ \psi' : \psi'^{-1}(U \cap U') \to \psi^{-1}(U \cap U')$, i.e. the change of variables map. Note that $y(x) = 0$, hence you have (a.e.) $$\partial_d(f \circ \psi')\big|_0 = \sum_{i=1}^d \partial_i(f \circ \psi)\big|_0 \partial_d y_i\big|_0.$$ Since it is certainly possible for the two sides of this equation to be different, it is meaningless to attempt to describe the vanishing of the normal derivative of $f$ in terms of the vanishing of the $d$-th partial derivative in generic boundary co-ordinates.

What is true is that, in each co-ordinate system $(U,\psi)$, the (unit) normal derivative of $f$ will correspond to a certain linear combination of partial derivatives of $f \circ \psi$. Notice that this linear combination will always feature a non-zero multiple of the $d$-th partial derivative.

If your boundary were better than Lipschitz (for instance if it were smooth) then you could also set up special coordinate systems adapted to $\partial D$ which, in addition to your requirements, also have the property that the $d$-th partial derivative of $f \circ \psi$ on $\{ x_d=0\}$ really does equal the (unit) normal derivative of $f$ relative to $\partial D$. These are known as geodesic (or Gaussian) normal co-ordinates (adapted to the hypersurface $\partial D$). I don't know what the theory is for these co-ordinate systems in such a low regularity setting as Lipschitz.

In order to answer your question you thus need some way of generating $C^2(\overline{D})$ functions $f$ with prescribed boundary data, i.e. with (a.e.) vanishing normal derivative together with the condition $f|_{\partial D} = g$ for some non-constant $g$. You should be able to do so using the celebrated Whitney extension theorem. I expect that the argument will go broadly speaking as follows: The vanishing of the normal derivative is essentially a constraint on a particular linear combination of partial derivatives of $f$. Assume WLOG that locally around $p \in \partial D$, the a.e. defined vector field of normal derivatives $\nu = (\nu_1, \ldots, \nu_d)$ is such that $\nu_d \neq 0$. Hence $$(\nu, \nabla f) = \left(\sum_{i=1}^{d-1} \nu_i \partial_i f\right) + \nu_d \partial_d f, \quad \text{and} \quad \partial_d f = -(\nu_d)^{-1}\sum_{i=1}^{d-1} \nu_i \partial_i f \ \text{if} \ (\nu, \nabla f) = 0.$$ Hence, you can choose anything you like for the partial derivatives $\partial_i f$, $i=1,\ldots,d-1$ on $\partial D$, and provided the $d$-th partial derivative is then given by the formula above, you have ensured that $f$ satisfies the Neumann boundary condition. Going the other way around, suppose you have the following collection of functions on $\partial D$: (i) a (non-constant, to solve your problem!) $f_0$ on $\partial D$; (ii) for $i=1, \ldots, d-1$, functions $f_i$; (iii) $f_d := -(\nu_d)^{-1} \sum_{i=1}^{d-1}\nu_i f_i$. Then (provided you can check that the hypotheses of the Whitney extension theorem are satisfied) there locally exists a function $f$, $C^2$ up to and including $\partial D$, which restricts to $f_0$ on $\partial D$ and which, by construction, satisfies the Neumann boundary condition there. From here, globalize to get the desired result.

This will only be a partial answer to your question, due to lack of time.

As you effectively noted, any $f\in C^2(\overline{D})$ with compact support and such that $(\nu, \nabla f) = 0 \ \mathcal{H}^{d-1}$-a.e. on $\partial D$ will belong to $D(L)$. The problem with your argument is that, given a differentiability point $p \in \partial D$ and a general $U \subseteq \mathbb{R}^d$ and $\psi : B(1) \to U$ as you describe, it is not true that $$(\nu, \nabla f) = 0 \ \text{at} \ p \Leftrightarrow \partial_d(f \circ \psi)\big|_0 = 0 .$$ To see why, suppose you had a different pair $(U',\psi')$ satisfying the same requirements. Then, by the chain rule, a.e. on $\psi'^{-1}(U \cap U')$ it holds that $$\partial_d(f \circ \psi')\big|_{x} = \partial_d(f \circ \psi \circ \psi^{-1} \circ \psi')\big|_{x} = \sum_{i=1}^d \partial_i(f \circ \psi)\big|_{y(x)}\partial_d y_i\big|_{x}$$ where $y := \psi^{-1} \circ \psi' : \psi'^{-1}(U \cap U') \to \psi^{-1}(U \cap U')$, i.e. the change of variables map. Note that $y(x) = 0$, hence you have (a.e.) $$\partial_d(f \circ \psi')\big|_0 = \sum_{i=1}^d \partial_i(f \circ \psi)\big|_0 \partial_d y_i\big|_0.$$ Since it is certainly possible for the two sides of this equation to be different, it is meaningless to attempt to describe the vanishing of the normal derivative of $f$ in terms of the vanishing of the $d$-th partial derivative in generic boundary co-ordinates.

What is true is that, in each co-ordinate system $(U,\psi)$, the (unit) normal derivative of $f$ will correspond to a certain linear combination of partial derivatives of $f \circ \psi$. Notice that this linear combination will always feature a non-zero multiple of the $d$-th partial derivative.

If your boundary were better than Lipschitz (for instance if it were smooth) then you could also set up special coordinate systems adapted to $\partial D$ which, in addition to your requirements, also have the property that the $d$-th partial derivative of $f \circ \psi$ on $\{ x_d=0\}$ really does equal the (unit) normal derivative of $f$ relative to $\partial D$. These are known as geodesic (or Gaussian) normal co-ordinates (adapted to the hypersurface $\partial D$). I don't know what the theory is for these co-ordinate systems in such a low regularity setting as Lipschitz.

In order to answer your question you thus need some way of generating $C^2(\overline{D})$ functions $f$ with prescribed boundary data, i.e. with (a.e.) vanishing normal derivative together with the condition $f|_{\partial D} = g$ for some non-constant $g$. You should be able to do so using the celebrated Whitney extension theorem. I expect that the argument will go broadly speaking as follows: The vanishing of the normal derivative is essentially a constraint on a particular linear combination of partial derivatives of $f$. Assume WLOG that locally around $p \in \partial D$, the a.e. defined vector field of normal derivatives $\nu = (\nu_1, \ldots, \nu_d)$ is such that $\nu_d \neq 0$. Hence $$(\nu, \nabla f) = \left(\sum_{i=1}^{d-1} \nu_i \partial_i f\right) + \nu_d \partial_d f, \quad \text{and} \quad \partial_d f = -(\nu_d)^{-1}\sum_{i=1}^{d-1} \nu_i \partial_i f \ \text{if} \ (\nu, \nabla f) = 0.$$ Hence, you can choose anything you like for the partial derivatives $\partial_i f$, $i=1,\ldots,d-1$ on $\partial D$, and provided the $d$-th partial derivative is then given by the formula above, you have ensured that $f$ satisfies the Neumann boundary condition. Going the other way around, suppose you have the following collection of functions on $\partial D$: (i) a (non-constant, to solve your problem!) $g$ on $\partial D$; (ii) for $i=1, \ldots, d-1$, functions $f_i$; (iii) $f_d := -(\nu_d)^{-1} \sum_{i=1}^{d-1}\nu_i f_i$. Then (provided you can check that the hypotheses of the Whitney extension theorem are satisfied) there locally exists a function $f$, $C^2$ up to and including $\partial D$, which restricts to $g$ on $\partial D$ and which, by construction, satisfies the Neumann boundary condition there. From here, globalize to get the desired result.

This will only be a partial answer to your question, due to lack of time.

As you effectively noted, any $f\in C^2(\overline{D})$ with compact support and such that $(\nu, \nabla f) = 0 \ \mathcal{H}^{d-1}$-a.e. on $\partial D$ will belong to $D(L)$. The problem with your argument is that, given a differentiability point $p \in \partial D$ and a general $U \subseteq \mathbb{R}^d$ and $\psi : B(1) \to U$ as you describe, it is not true that $$(\nu, \nabla f) = 0 \ \text{at} \ p \Leftrightarrow \partial_d(f \circ \psi)\big|_0 = 0 .$$ To see why, suppose you had a different pair $(U',\psi')$ satisfying the same requirements. Then, by the chain rule, a.e. on $\psi'^{-1}(U \cap U')$ it holds that $$\partial_d(f \circ \psi')\big|_{x} = \partial_d(f \circ \psi \circ \psi^{-1} \circ \psi')\big|_{x} = \sum_{i=1}^d \partial_i(f \circ \psi)\big|_{y(x)}\partial_d y_i\big|_{x}$$ where $y := \psi^{-1} \circ \psi' : \psi'^{-1}(U \cap U') \to \psi^{-1}(U \cap U')$, i.e. the change of variables map. Note that $y(x) = 0$, hence you have (a.e.) $$\partial_d(f \circ \psi')\big|_0 = \sum_{i=1}^d \partial_i(f \circ \psi)\big|_0 \partial_d y_i\big|_0.$$ Since it is certainly possible for the two sides of this equation to be different, it is meaningless to attempt to describe the vanishing of the normal derivative of $f$ in terms of the vanishing of the $d$-th partial derivative in generic boundary co-ordinates.

What is true is that, in each co-ordinate system $(U,\psi)$, the (unit) normal derivative of $f$ will correspond to a certain linear combination of partial derivatives of $f \circ \psi$. Notice that this linear combination will always feature a non-zero multiple of the $d$-th partial derivative.

If your boundary were better than Lipschitz (for instance if it were smooth) then you could also set up special coordinate systems adapted to $\partial D$ which, in addition to your requirements, also have the property that the $d$-th partial derivative of $f \circ \psi$ on $\{ x_d=0\}$ really does equal the (unit) normal derivative of $f$ relative to $\partial D$. These are known as geodesic (or Gaussian) normal co-ordinates (adapted to the hypersurface $\partial D$). I don't know what the theory is for these co-ordinate systems in such a low regularity setting as Lipschitz.

In order to answer your question you thus need some way of generating $C^2(\overline{D})$ functions $f$ with prescribed boundary data, i.e. with (a.e.) vanishing normal derivative together with the condition $f|_{\partial D} = g$ for some non-constant $g$. You should be able to do so using the celebrated Whitney extension theorem. I expect that the argument will go broadly speaking as follows: The vanishing of the normal derivative is essentially a constraint on a particular linear combination of partial derivatives of $f$. Assume WLOG that locally around $p \in \partial D$, the a.e. defined vector field of normal derivatives $\nu = (\nu_1, \ldots, \nu_d)$ is such that $\nu_d \neq 0$. Hence $$(\nu, \nabla f) = \left(\sum_{i=1}^{d-1} \nu_i \partial_i f\right) + \nu_d \partial_d f, \quad \text{and} \quad \partial_d f = -(\nu_d)^{-1}\sum_{i=1}^{d-1} \nu_i \partial_i f \ \text{if} \ (\nu, \nabla f) = 0.$$ Hence, you can choose anything you like for the partial derivatives $\partial_i f$, $i=1,\ldots,d-1$ on $\partial D$, and provided the $d$-th partial derivative is then given by the formula above, you have ensured that $f$ satisfies the Neumann boundary condition. Going the other way around, suppose you have the following collection of functions on $\partial D$: (i) a (non-constant, to solve your problem!) $f_0$ on $\partial D$; (ii) for $i=1, \ldots, d-1$, functions $f_i$; (iii) $f_d := -(\nu-d)^{-1} \sum_{i=1}^{d-1}\nu_i f_i$. Then (provided you can check that the hypotheses of the Whitney extension theorem are satisfied) there locally exists a function $f$, $C^2$ up to and including $\partial D$, which restricts to $f_0$ on $\partial D$ and which, by construction, satisfies the Neumann boundary condition there. From here, globalize to get the desired result.

This will only be a partial answer to your question, due to lack of time.

As you effectively noted, any $f\in C^2(\overline{D})$ with compact support and such that $(\nu, \nabla f) = 0 \ \mathcal{H}^{d-1}$-a.e. on $\partial D$ will belong to $D(L)$. The problem with your argument is that, given a differentiability point $p \in \partial D$ and a general $U \subseteq \mathbb{R}^d$ and $\psi : B(1) \to U$ as you describe, it is not true that $$(\nu, \nabla f) = 0 \ \text{at} \ p \Leftrightarrow \partial_d(f \circ \psi)\big|_0 = 0 .$$ To see why, suppose you had a different pair $(U',\psi')$ satisfying the same requirements. Then, by the chain rule, a.e. on $\psi'^{-1}(U \cap U')$ it holds that $$\partial_d(f \circ \psi')\big|_{x} = \partial_d(f \circ \psi \circ \psi^{-1} \circ \psi')\big|_{x} = \sum_{i=1}^d \partial_i(f \circ \psi)\big|_{y(x)}\partial_d y_i\big|_{x}$$ where $y := \psi^{-1} \circ \psi' : \psi'^{-1}(U \cap U') \to \psi^{-1}(U \cap U')$, i.e. the change of variables map. Note that $y(x) = 0$, hence you have (a.e.) $$\partial_d(f \circ \psi')\big|_0 = \sum_{i=1}^d \partial_i(f \circ \psi)\big|_0 \partial_d y_i\big|_0.$$ Since it is certainly possible for the two sides of this equation to be different, it is meaningless to attempt to describe the vanishing of the normal derivative of $f$ in terms of the vanishing of the $d$-th partial derivative in generic boundary co-ordinates.

What is true is that, in each co-ordinate system $(U,\psi)$, the (unit) normal derivative of $f$ will correspond to a certain linear combination of partial derivatives of $f \circ \psi$. Notice that this linear combination will always feature a non-zero multiple of the $d$-th partial derivative.

If your boundary were better than Lipschitz (for instance if it were smooth) then you could also set up special coordinate systems adapted to $\partial D$ which, in addition to your requirements, also have the property that the $d$-th partial derivative of $f \circ \psi$ on $\{ x_d=0\}$ really does equal the (unit) normal derivative of $f$ relative to $\partial D$. These are known as geodesic (or Gaussian) normal co-ordinates (adapted to the hypersurface $\partial D$). I don't know what the theory is for these co-ordinate systems in such a low regularity setting as Lipschitz.

In order to answer your question you thus need some way of generating $C^2(\overline{D})$ functions $f$ with prescribed boundary data, i.e. with (a.e.) vanishing normal derivative together with the condition $f|_{\partial D} = g$ for some non-constant $g$. You should be able to do so using the celebrated Whitney extension theorem. I expect that the argument will go broadly speaking as follows: The vanishing of the normal derivative is essentially a constraint on a particular linear combination of partial derivatives of $f$. Assume WLOG that locally around $p \in \partial D$, the a.e. defined vector field of normal derivatives $\nu = (\nu_1, \ldots, \nu_d)$ is such that $\nu_d \neq 0$. Hence $$(\nu, \nabla f) = \left(\sum_{i=1}^{d-1} \nu_i \partial_i f\right) + \nu_d \partial_d f, \quad \text{and} \quad \partial_d f = -(\nu_d)^{-1}\sum_{i=1}^{d-1} \nu_i \partial_i f \ \text{if} \ (\nu, \nabla f) = 0.$$ Hence, you can choose anything you like for the partial derivatives $\partial_i f$, $i=1,\ldots,d-1$ on $\partial D$, and provided the $d$-th partial derivative is then given by the formula above, you have ensured that $f$ satisfies the Neumann boundary condition. Going the other way around, suppose you have the following collection of functions on $\partial D$: (i) a (non-constant, to solve your problem!) $f_0$ on $\partial D$; (ii) for $i=1, \ldots, d-1$, functions $f_i$; (iii) $f_d := -(\nu_d)^{-1} \sum_{i=1}^{d-1}\nu_i f_i$. Then (provided you can check that the hypotheses of the Whitney extension theorem are satisfied) there locally exists a function $f$, $C^2$ up to and including $\partial D$, which restricts to $f_0$ on $\partial D$ and which, by construction, satisfies the Neumann boundary condition there. From here, globalize to get the desired result.