# Reproducing figure from paper on hydrodynamics

I am trying to solve Eqn. on Page 674 of https://www.nature.com/articles/nphys3667.pdf I shall quote the relevant paragraph :

Next, we proceed to analyse the nonlocal response in a strip with transverse current injected and drained through a pair of contacts, as pictured in Fig. 1(a). Unlike the above case of longitudinal current, here the potential profile is not set externally but must be obtained from equation (3). The analysis is facilitated by introducing a stream function via $$v = \hat{z}$$ × ∇ $$\psi$$ , which solves the incompressibility condition. At first we will completely ignore the ohmic effects, setting α and γp to zero as above, which leads to a biharmonic equation $$(\partial_x^2+\partial_y^2)^2\psi = 0$$ with the boundary conditions $$v_x$$ = 0, $$v_y$$ =$$I$$δ(x) for y =$$0,w$$.

The geometry is $$(x,y)\in [0,l]\times[0,w]$$ and $$l \gg w$$. The equation to solve is $$\nabla^2\nabla^2\psi = 0$$ where $$\psi$$ is the stream function. I write this as two coupled laplacians, $$\nabla^2\psi = u$$ and $$\nabla^2 u = 0$$. The boundary conditions are current injected at $$y=0,w$$. The current injected is $$I = n e \delta(x)$$, and I can model the delta as a lorentzian, so $$I = n e v_y = n e \frac{1}{(x/a)^2+1}$$, i.e., $$v_y = \frac{1}{(x/a)^2+1}$$, and since $$v = \hat{z}\times\partial\psi = (-\partial_y\psi, \partial_x\psi)$$, we have $$\partial_x\psi = \frac{1}{(x/a)^2+1}$$, thus $$\psi = \tan^{-1}(x/a)$$ at $$y=0,w$$ and $$u = \nabla^2\psi = \frac{-2x}{a^2}\frac{1}{(x/a)^2+1}$$ at $$y = 0,w$$. At $$x=0,l$$ I have set $$\psi$$ constant at $$x=0,l$$.

The physical picture that one should get is Figure 1 of the paper, which has vortices, which looks like : I cannot seem to reproduce these vortices. I model the $$\delta(x) = \frac{1}{x^2+1}$$. My code is below :

L1 = 35 ; H1 = 5  ; \[CapitalOmega]1 =
Rectangle[{-L1/2, -H1/2}, {L1/2, H1/
2}] ; RegionPlot[\[CapitalOmega]1, AspectRatio -> Automatic]

mesh1 = ToElementMesh[ \[CapitalOmega]1, MaxCellMeasure -> 0.001];

biharmonic = { Laplacian[psi2[x, y], {x, y}] == u2[x, y],
Laplacian[u2[x, y], {x, y}] == 0};

bcond = {DirichletCondition[psi2[x, y] == ArcTan[x], y == -H1/2],
DirichletCondition[psi2[x, y] ==  ArcTan[x], y == H1/2],
DirichletCondition[u2[x, y] == -((2 x)/(x^2 + 1)^2), y == -H1/2],
DirichletCondition[u2[x, y] == -((2 x)/(x^2 + 1)^2) , y == H1/2],
DirichletCondition[psi2[x, y] ==  1, x == -L1/2],
DirichletCondition[psi2[x, y] ==  -1, x == L1/2]};

Needs["NDSolveFEM"]
soln = NDSolve[{biharmonic, bcond}, {psi2, u2},
Element[{x, y}, mesh1]];

StreamDensityPlot[
Evaluate[{-Derivative[0, 1][psi2][x, y],
Derivative[1, 0][psi2][x, y]} /. soln[]],
Element[{x, y}, mesh1], PlotLegends -> Automatic,
AspectRatio -> Automatic, ColorFunction -> cool, ImageSize -> Large]


Experts, am I doing something wrong ? Thanks.

• Include Needs["NDSolveFEM"] before mesh1 = ToElementMesh... May 12 at 7:34
• What were the error messeges that you got when you ran the code?
– YSP
May 12 at 8:29
• @YSP No error messages, does not match physical picture May 12 at 8:33
• @DanielHuber Thanks. Good catch, I was copy pasting blocks of code and well...but I don't think that was the problem because I had already loaded the package on this notebook. May 12 at 8:35

We can't reproduce picture from the paper cited since they used Fourie integral on $$-\infty \le x \le \infty$$. With FEM we can reproduce part of picture be crossing side vortexes as follows

Clear["Global*"]
Needs["NDSolveFEM"]

L1 = 2; H1 = 1; \[CapitalOmega]1 =
Rectangle[{-L1/2, -H1/2}, {L1/2, H1/2}];

mesh1 = ToElementMesh[\[CapitalOmega]1,
MaxCellMeasure -> 0.002]; mesh1["Wireframe"]

biharmonic = {Laplacian[psi2[x, y], {x, y}] == u2[x, y],
Laplacian[u2[x, y], {x, y}] == 0};

bcond = {DirichletCondition[psi2[x, y] == -UnitStep[x], y == -H1/2],
DirichletCondition[psi2[x, y] == -UnitStep[x], y == H1/2],
DirichletCondition[{psi2[x, y] == -y^2 + (H1/2)^2, u2[x, y] == 0},
x == -L1/2],
DirichletCondition[{psi2[x, y] == y^2 - (H1/2)^2 - 1,
u2[x, y] == 0}, x == L1/2]};

soln = NDSolve[{biharmonic, bcond}, {psi2, u2},
Element[{x, y}, mesh1]];

Show[DensityPlot[
Evaluate[Norm[{-Derivative[0, 1][psi2][x, y],
Derivative[1, 0][psi2][x, y]} /. soln[]]],
Element[{x, y}, mesh1], PlotLegends -> Automatic,
ColorFunction -> "MintColors", PlotPoints -> 50,
AspectRatio -> Automatic, Frame -> False, PlotRange -> {0, 2}],
StreamPlot[
Evaluate[{-Derivative[0, 1][psi2][x, y],
Derivative[1, 0][psi2][x, y]} /. soln[]],
Element[{x, y}, mesh1], StreamPoints -> Fine,
StreamColorFunction -> None, StreamStyle -> Magenta]] Update 1. We also can reproduce vortexes by adding two logarithmic singularities as follows

Clear["Global*"]
Needs["NDSolveFEM"]

L1 = 3; H1 = 1; \[CapitalOmega]1 =
Rectangle[{-L1/2, -H1/2}, {L1/2, H1/2}];

mesh1 = ToElementMesh[\[CapitalOmega]1,
MaxCellMeasure -> 0.002]; mesh1["Wireframe"]

biharmonic = {Laplacian[psi2[x, y], {x, y}] == u2[x, y],
Laplacian[u2[x, y], {x, y}] == 0};

bcond = {DirichletCondition[psi2[x, y] == -UnitStep[x], y == -H1/2],
DirichletCondition[psi2[x, y] == -UnitStep[x], y == H1/2],
DirichletCondition[{psi2[x, y] == -Log[(x + H1)^2 + y^2],
u2[x, y] == 0}, x == -L1/2],
DirichletCondition[{psi2[x, y] == Log[(x - H1)^2 + y^2] - 1,
u2[x, y] == 0}, x == L1/2]};

soln = NDSolve[{biharmonic, bcond}, {psi2, u2},
Element[{x, y}, mesh1]];

Show[DensityPlot[
Evaluate[Norm[{-Derivative[0, 1][psi2][x, y],
Derivative[1, 0][psi2][x, y]} /. soln[]]],
Element[{x, y}, mesh1], PlotLegends -> Automatic,
ColorFunction -> "MintColors", PlotPoints -> 50,
AspectRatio -> Automatic, Frame -> False, PlotRange -> {0, 2}],
StreamPlot[
Evaluate[{-Derivative[0, 1][psi2][x, y],
Derivative[1, 0][psi2][x, y]} /. soln[]],
Element[{x, y}, mesh1], StreamPoints -> Fine,
StreamColorFunction -> None, StreamStyle -> Magenta]] Update 2. We also can use Fourie series as follows

 Clear[p]

nmax = 50; L = 65;
psi[x_, y_] := Sum[Sin[ Pi n/L  x] p[n][y], {n, nmax}];
k[n_] := Pi n/L;
psix[x_, y_] := Sum[k[n] Cos[k[n]  x] p[n][y], {n, nmax}];

c[n_] = FourierSinCoefficient[-UnitStep[x], x, n]

sol1 = DSolve[{k[n]^4 - 2 k[n]^2 D[pn[y], {y, 2}] + D[pn[y], {y, 4}] ==
0, pn[-1/2] == c[n], pn[1/2] == c[n], pn'[-1/2] == 0,
pn'[1/2] == 0}, pn[y], y] // FullSimplify;

p[m_][y_] := (pn[y] /. sol1[]) /. n -> m

p1[n_][y_] := D[pn[y] /. sol1[], y]

psiy[x_, y_] := Sum[Sin[k[n]  x] p1[n][y], {n, nmax}];

StreamDensityPlot[
Evaluate[{-psiy[x, y], psix[x, y]}], {x, -1.5, 1.5}, {y, -.5, .5},
AspectRatio -> Automatic, StreamPoints -> Fine,
ColorFunction -> "MintColors"] • Alex, that makes sense. Guess need to do something like Show[MapAt[Translate[#, {-L1, 0}] &, plot, 1], plot, MapAt[Translate[#, {L1, 0}] &, plot, 1], PlotRange -> All] to get a periodic picture. Let me think about this a bit more before I accept your answer, thanks. May 12 at 9:49
• @Charlie Please, see Update 2 to my answer. May 13 at 14:42