# Finite-Slip Boundary Conditions

## Equations

I am trying to solve the linearized Navier-Stokes equations for a 2D electron fluid given by:

\begin{aligned} \sigma \nabla \Phi + D^2 \nabla^2 \boldsymbol{J}&=\boldsymbol{J} \\ \nabla \cdot \boldsymbol{J} &=0, \end{aligned}

where $$\sigma$$ is the conductivity, and $$D$$ is a parameter called the Gurzhi length which is related to the diffusion-length of vorticity in the system (Note that as $$D\rightarrow0$$, we recover Ohm's law).

First, let's express this in the time-independent coefficients form:

$$\nabla \cdot \left(-c\nabla \boldsymbol{u} - \alpha \boldsymbol{u} + \gamma \right) + \beta \cdot \nabla \boldsymbol{u} + a \boldsymbol{u} = f,$$

where we've identified the field and coefficients as:

$$\boldsymbol{u}=\begin{pmatrix} \Phi \\ J_x \\ J_y \end{pmatrix}$$, $$a = \left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right)$$, $$c=\left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & d^2 & 0 \\ 0 & 0 & d^2 \\ \end{array} \right)$$, $$\beta = \left( \begin{array}{ccc} \begin{pmatrix}0 \\0\end{pmatrix} & \begin{pmatrix}1 \\0\end{pmatrix} & \begin{pmatrix}0 \\1\end{pmatrix} \\ \begin{pmatrix}\sigma \\0\end{pmatrix} & \begin{pmatrix}0 \\0\end{pmatrix} & \begin{pmatrix}0 \\0\end{pmatrix} \\ \begin{pmatrix}0 \\\sigma\end{pmatrix} & \begin{pmatrix}0 \\0\end{pmatrix} & \begin{pmatrix}0 \\0\end{pmatrix} \\ \end{array} \right)$$, $$\alpha=\gamma=f=0$$

and use the PDE building blocks introduced in version 13 for convenience:

vars = {ϕ[x, y], Jx[x, y], Jy[x, y]};
reactionTerm = ReactionPDETerm[{vars, {x, y}},
DiagonalMatrix[{0, 1, 1}]];

diffusionTerm[d_] = DiffusionPDETerm[{vars, {x, y}},
{{0, 0, 0}, {0, {d^2, d^2}, 0}, {0, 0, {d^2, d^2}}}];

convectionTerm[σ_] = ConvectionPDETerm[{vars, {x, y}},
{{{0, 0}, {1, 0}, {0, 1}},
{{σ, 0}, {0, 0}, {0, 0}},
{{0, σ}, {0, 0}, {0, 0}}}];

partialDifferentialEquation[d_, σ_] = reactionTerm + diffusionTerm[d] + convectionTerm[σ]


## No-Slip Boundary Conditions

We can solve this using appropriate boundary conditions. E.g. consider DirichletConditions for $$\Phi$$, and the common no-slip boundary conditions for $$J_x,J_y$$, in a simple channel flow geometry:

reg = Rectangle[{-1/2, -1/2}, {1/2, 1/2}];
noSlipBCs = {
DirichletCondition[ϕ[x, y] == 0., x == 1/2],
DirichletCondition[ϕ[x, y] == 1., x == -1/2],
DirichletCondition[{Jx[x, y] == 0., Jy[x, y] == 0.},
y == -1/2 || y == 1/2]
};

parametricSolutions = ParametricNDSolveValue[{
partialDifferentialEquation[d, 1] == {0, 0, 0},
noSlipBCs}, {ϕ, Jx, Jy}, {x, y} ∈ reg, {d},
Method -> {"FiniteElement",
"InterpolationOrder" -> {ϕ -> 1, Jx -> 2, Jy -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}]


Varying $$D$$ and plotting the normalized horizontal current density, we indeed recover something resembling uniform current density as $$D\rightarrow 0$$, and parabolic current density as $$D\rightarrow \infty$$:

normalizedHorizontalCurrent[d_] :=
Module[{sol = parametricSolutions[d][[2]], totalCurrent},
totalCurrent = NIntegrate[sol[0., y], {y, -1/2, 1/2}];
sol[0., #]/totalCurrent &]

With[{sols = normalizedHorizontalCurrent /@ {0.01, 0.1, 1}},
Plot[Through[sols[y]] // Evaluate, {y, -1/2, 1/2}, Frame -> True,
PlotRange -> {0, 3/2}, PlotLegends -> {0.01, 0.1, 1}]]


## Finite-Slip Boundary Conditions?

We now wish to impose finite-slip boundary conditions as described in equations (1-3) here, given by:

$$\boldsymbol{J}^t = \zeta \hat{n} \cdot \nabla \boldsymbol{J}^t,$$

where $$\boldsymbol{J}^t= \boldsymbol{J} - \left(\boldsymbol{J}\cdot \hat{n}\right)\hat{n}$$ is the tangential current density at the boundary, and $$\zeta$$ is the slip length. Note that as $$\zeta \rightarrow 0$$ and $$\zeta \rightarrow \infty$$ we recover the no-slip and no-stress boundary conditions respectively.

The question is: how do we specify this boundary condition using NeumannValue? A comment by user21 in this related question suggests we can use a combination of Indexed and BoundaryUnitNormal to specify this.

Note: I'm aware the simple geometry has both a closed form solution, and that the boundary normals are particularly simple. The actual geometries I want to study do not, hence a general solution in terms of BoundaryUnitNormal (or similar) is highly desired.

### EDIT 01

First, simplify expression using {nx,ny}:

With[{n = {nx, ny}},
FullSimplify[
n .(Grad[{Jx[x, y], Jy[x, y]} - ({Jx[x, y], Jy[x, y]} . n) n, {x,y}])]
]


{-((-1+nx^2+ny^2) (nx (Jx^(1,0))[x,y]+ny (Jy^(1,0))[x,y])),-((-1+nx^2+ny^2) (nx (Jx^(0,1))[x,y]+ny (Jy^(0,1))[x,y]))}

And then replace {nx,ny} with BoundaryUnitNormal:

Needs["NDSolveFEM"]
neumann = With[{n = BoundaryUnitNormal[x, y]},
-(n . n - 1) n . Grad[{Jx[x, y], Jy[x, y]}, {x, y}]]

finiteSlipNeumann = {0,
NeumannValue[Indexed[0.1 neumann, 1], y == -1/2 || y == 1/2],
NeumannValue[Indexed[0.1 neumann, 2], y == -1/2 || y == 1/2]
};


By removing the DirichletCondition on Jx and Jy, and adding finiteSlipNeumann to the RHS of the equation inside ParametricNDSolveValue, but I get this error:

Derivatives of dependent variables in boundary conditions are not supported with the Finite Element Method in this version of NDSolve.

What is more, the expression I derived above is certainly wrong, since $$\hat{n}$$ is normalized and thus $$\hat{n}\cdot\hat{n}-1$$ evaluates to zero.

• I have some trouble understanding your equation setup. If sigma is conductivity, should this then not go in a diffusive therm (c)? If we assume that the equation is correct (\sigma \nabla \phi), then \sigma needs to be a vector (sigmaX, sigmaY) for it to multiply with gradient phi. Also I do not see terms grad Jx and grad Jy in your equations, yet you model them. Commented Feb 14, 2022 at 9:26
• Concerning BoundaryUnitNormal this represents {nx, ny} and with Indexed you can extract a part of that vector, if needed. Commented Feb 14, 2022 at 9:28
• Thanks @user21 - The equation is correct (and indeed \sigma here is isotropic), although I forgot to specify the continuity equation (essentially J acts as an incompressible fluid). See edit. Re: the NeumannValue - I'm not sure I expressed it correctly, since I'm getting an error about derivatives of dependent variables.. Commented Feb 14, 2022 at 13:54
• In-fact, pde error aside - the expression I used is certainly wrong since the normals are normalized and thus {nx,ny}.{nx,ny}-1 evaluates to zero Commented Feb 14, 2022 at 14:22
• Currently (V13.0), you can not use the derivative of the dependent variable in NeumannValue. But this came up the third time now the past week. Seems relevant to have.... Commented Feb 14, 2022 at 15:46