4
$\begingroup$

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}]]

enter image description here

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["NDSolve`FEM`"]
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.

$\endgroup$
5
  • 1
    $\begingroup$ 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. $\endgroup$
    – user21
    Commented Feb 14, 2022 at 9:26
  • 1
    $\begingroup$ Concerning BoundaryUnitNormal this represents {nx, ny} and with Indexed you can extract a part of that vector, if needed. $\endgroup$
    – user21
    Commented Feb 14, 2022 at 9:28
  • $\begingroup$ 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.. $\endgroup$ Commented Feb 14, 2022 at 13:54
  • $\begingroup$ 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 $\endgroup$ Commented Feb 14, 2022 at 14:22
  • 1
    $\begingroup$ 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.... $\endgroup$
    – user21
    Commented Feb 14, 2022 at 15:46

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.