# Why does NDSolve solution not match physical intuition?

Related to a previous question I posted, I am trying to solve for the electric potential on a box which obeys Ohm's law: $$$$\nabla\cdot(\overset{\scriptscriptstyle\leftrightarrow}{\sigma} \nabla \Phi) = 0, \label{Eqn:OhmsLaw}$$$$ where $$\overset{\scriptscriptstyle\leftrightarrow}{\sigma}$$ is the 2D conductivity tensor, and $$\Phi$$ is the electric potential. The diagonal and off-diagonal elements of the conductivity tensor are given by the functions defined in the Mathematica code.

When $$\sigma_{xx}$$ is zero in an interval Ohms law reduces to $$\hat{z}\cdot((\nabla\sigma_{xy})\times(\nabla\Phi))$$, thus contours of electric potential are parallel to contour of constant $$\sigma_{xy}$$. Since $$\sigma_{xy}$$ is radially symmetric equipotential curves of $$\Phi$$ should be circles.

Here is the Mathematica code :

ClearAll["Global*"]
<< NDSolveFEM

nl = 40 ;
n1 = 2 ;
rpn = 0.6 ;
delta = 0.01 ;
w = 0.5 ;
r0 = rpn + delta + w ;
r1 = rpn - delta  - w ;
a = Solve[n1 == nl (1 - a/r0^3), a][[1]][[1]][[2]];
N1 = 40 (1 - a/r^3);
N2 = 2;
N3 = 2/delta (r - rpn) ;
N4  = -2 ;
beta1 = Solve[1 - 1/beta == 5 (1 - 1/(r1^2 + beta^2)^(1/2)),
beta ][[3]][[1]][[2]];
gamma = -10/(1 - 1/beta1);
N5 = gamma (1 - 1/Sqrt[r^2 + beta1^2]);
micrometer = 1;
centimeter = micrometer 10^4;
gauss = 1 ;
tesla = 10^4 gauss ;
Lx = 15; Ly = 15;
omegactimestau[B_] := 4.8 10^-10 B/(1.7 10^-29 3 10^10) 3 10^-4/10^8;

\[Sigma]xx = Piecewise[{
{Abs[N5 + 2], 0 <= r <= r1  },
{0 , r1 <= r <= rpn - delta},
{Sqrt[(2 + N3)] Sqrt[(2 - N3)], rpn - delta <= r <= rpn + delta},
{0 , rpn + delta <= r <= r0},
{N1 - 2, r >= r0}
}];
\[Sigma]xy =  Piecewise[{
{N5, 0 <= r <= r1  },
{N4, r1 <= r <= rpn - delta},
{N3, rpn - delta <= r <= rpn + delta},
{N2, rpn + delta <= r <= r0},
{N1, r >= r0}
}];;

sigmaxxcheck[x_, y_, B_] :=
1/omegactimestau[B] \[Sigma]xx /. {r ->
Sqrt[(x - Lx/2)^2 + (y - Ly/2)^2]}
sigmaxycheck[x_, y_,
B_] := \[Sigma]xy /. {r -> Sqrt[(x - Lx/2)^2 + (y - Ly/2)^2]}
rectangle = ImplicitRegion[0 <= x <= Lx && 0 <= y <= Ly, {x, y}];
frefine =
Function[{vertices, area}, Block[{x1, y1}, {x1, y1} = Mean[vertices];
If[(x1 - Lx/2)^2 + (y1 - Ly/2)^2 < 4 && area > 0.0005, True,
False]]];

mesh = ToElementMesh[rectangle, MaxBoundaryCellMeasure -> 10^-3,
MeshRefinementFunction -> frefine];
reg = mesh;

brange = Range[0.5, 2, 0.5]*tesla;
sol = NDSolveValue[{Inactive[
Div][{{sigmaxxcheck[x, y, #],
sigmaxycheck[x, y, #]}, {-sigmaxycheck[x, y, #],
sigmaxxcheck[x, y, #]}} .
Inactive[Grad][phi[x, y], {x, y}], {x, y}] == 0,
phi[Lx, y] == 10, phi[0, y] == -10},
phi, {x, y} \[Element] reg] & /@ brange;

ContourPlot[
sol[x + Lx/2, y + Ly/2] /. sol -> sol[[2]], {x, -1, 1}, {y, -1, 1},
AspectRatio -> Automatic, ColorFunction -> ColorData["SunsetColors"],
ImageSize -> 650, PlotRangePadding -> None, FrameStyle -> Black,
BaseStyle -> FontSize -> 22, PerformanceGoal -> "Quality",
Contours -> 100, PlotRange -> All, Exclusions -> None]
`

What I obtain is a mess.

Does this have something to do with the system being hyperbolic in some region and elliptic in another?

• Can include the code for the plot and the plot itself in the post? Jun 29, 2023 at 3:25
• @user21 I have included the code for the plot Jun 29, 2023 at 3:43
• And how about the plot? It's just a copy from Mathematica and paste in SE Jun 29, 2023 at 4:04
• @user21 I have edited to include plot output from Mathematica. Jun 29, 2023 at 4:18
• When I copy and paste this code I get error messages. Have you tried to load the FEM package right at the beginning. You'd need to call the package before you make use of it. Jun 29, 2023 at 4:45

With the explicit differential equation missing, only generalities can be cited (Citation by whome?: Every blockhead can cite generalities, it is the mastermind who differentiates all the special cases they constitute)

The physical intuition is another one and has nothing to do with symmetries.

$$\nabla \cdot \sigma \cdot \vec v == \nabla \cdot (\sigma \cdot \nabla \ \Phi) =0$$

is nothing more, than the conservation law of charge and plus a rotation free current having a (scalar) potential.

Your $$\nabla\times\nabla$$ from a force orthogonal to the gradient realizes
$$\nabla\times\nabla\times \vec v = \nabla . (\nabla \cdot \vec v ) -\Delta \vec v$$ the Laplacian for vector fields.

The vector Laplacian is, locally, defined as the volume mean integral, defined by Stokes theorem as the limit of the outward flow over the surface.

Flow orthogonal to lines of constant potential is common and known from wind on the rotating earth (Coriolis force) or 2d-electric current in magnetic fields (Hall effect).

Except for (mass or charge) conservation, there is no more information without boundary conditons and start field confiruration for time dependence.

In 2d, the equations $$\nabla \cdot \vec v =0$$ and $$\nabla \times \vec v ==$$ are identical with Cauchy-Riemann equations for the complex differential $$\partial_{z^\star}f(z) =0$$

So at the end, there is perhaps an obstacle: To describe a rotation free current as the gradient of a scalar potential. In electrodynamics, generally, all kinds of vector fields are components of exterior derivatives of other vector fields, the electric potential e.g. is the time component of the vector potential.

• How exactly is this useful? Jun 29, 2023 at 20:07
• In the same way as the question. By 99/100 Mathematica problems of physical origin arise from the difficult mathematical background problems implemented without much knowlegde of the theory. To answer your question I need to hold a one term course in Navier Stokes equations of plasma physics. Charlie will understand. Jun 30, 2023 at 4:44
• Oh interesting, how is the Navier-Stokes equation, or plasma physics involved? We should try a magneto-hydro-dymaics (MHD) approach! Jun 30, 2023 at 9:56