System of pde with Neumann boundary conditions

Question

The Ginzburg-Landau equation for a system of squared superconductor with a slit in presence of a uniform magnetic field can be written as a sort of continuity equation with Neumann boundary conditions (Ref):

$d_a. \frac{\partial \bf{u}}{\partial t} + \nabla. \Gamma = \bf{F}$

$\nabla \begin{pmatrix} u_1 \\ u_2 \end{pmatrix}.\bf{n}= 0$

$\nabla \times \begin{pmatrix} u_3 \\ u_4 \end{pmatrix} = \bf{B_a}$

$\begin{pmatrix} u_3 \\ u_4 \end{pmatrix}. \bf{n} = 0$

The magnetic field is applied along the z-axis so $\bf{B_a} = B_a \hat{k}, \bf{n} = \hat{k}$, and $B_a$ is constant. The result is a set of four coupled pde :

dA = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, sigma, 0}, {0, 0, 0, sigma}};

u[t_?NumericQ, x_?NumericQ, y_?NumericQ, z_?NumericQ] := {u1[t, x, y, z], u2[t, x, y, z], u3[t, x, y, z], 
  u4[t, x, y, z]}

delgma[t_, x_, y_, z_] = {(-D[u1[t, x, y, z], {x, 2}] - 
      D[u1[t, x, y, z], {y, 2}])/(kappa^2), (-D[
        u2[t, x, y, z], {x, 2}] - 
      D[u2[t, x, y, z], {y, 2}])/(kappa^2), 
   D[D[u4[t, x, y, z], x], y] - D[u3[t, x, y, z], {y, 2}] + 
    D[ba[x, y, z], y], 
   D[D[u3[t, x, y, z], x], y] - D[u4[t, x, y, z], {x, 2}] + 
    D[ba[x, y, z], x]};

f1[t_, x_, y_, z_] = -(D[u3[t, x, y, z], x] + D[u4[t, x, y, z], y])*
    u1[t, x, y, z]/
     kappa - (D[u1[t, x, y, z], x]*u3[t, x, y, z] + 
      D[u1[t, x, y, z], y]*u4[t, x, y, z])/
    kappa - (u3[t, x, y, z]*u3[t, x, y, z] + 
      u4[t, x, y, z]*u4[t, x, y, z])*u1[t, x, y, z] + (1 - temp)*
    u1[t, x, y, 
     z]*(1 - (u1[t, x, y, z]*u1[t, x, y, z] + 
        u2[t, x, y, z]*u2[t, x, y, z]));

f2[t_, x_, y_, z_] = -(D[u3[t, x, y, z], x] + D[u4[t, x, y, z], y])*
    u2[t, x, y, z]/
     kappa - (D[u2[t, x, y, z], x]*u3[t, x, y, z] + 
      D[u2[t, x, y, z], y]*u4[t, x, y, z])/
    kappa - (u3[t, x, y, z]*u3[t, x, y, z] + 
      u4[t, x, y, z]*u4[t, x, y, z])*u2[t, x, y, z] + (1 - temp)*
    u2[t, x, y, 
     z]*(1 - (u1[t, x, y, z]*u1[t, x, y, z] + 
        u2[t, x, y, z]*u2[t, x, y, z]));

f3[t_, x_, y_, z_] = (D[u2[t, x, y, z], x]*u1[t, x, y, z] - 
      D[u1[t, x, y, z], x]*u2[t, x, y, z])/
    kappa - (u1[t, x, y, z]*u1[t, x, y, z] + 
      u2[t, x, y, z]*u2[t, x, y, z])*u3[t, x, y, z];

f4[t_, x_, y_, z_] = (D[u2[t, x, y, z], y]*u1[t, x, y, z] - 
      D[u1[t, x, y, z], y]*u2[t, x, y, z])/
    kappa - (u1[t, x, y, z]*u1[t, x, y, z] + 
      u2[t, x, y, z]*u2[t, x, y, z])*u3[t, x, y, z];

f[t_, x_, y_, z_] = {f1[t, x, y, z], f2[t, x, y, z], f3[t, x, y, z], 
   f4[t, x, y, z]};

dudt[t_, x_, y_, z_] = {D[u1[t, x, y, z], t], D[u2[t, x, y, z], t], 
   D[u3[t, x, y, z], t], D[u4[t, x, y, z], t]};

eqns = dA.dudt[t, x, y, z] + delgma[t, x, y, z] - f[t, x, y, z];

The equations have to be solved over a squared region with a slit which I have tried to model as:

region3d = 
  ImplicitRegion[-1 <= x <= 1 && -1 <= y <= 
     1 && ! (0 < x <= 1 && -0.1 < y < 0.1) && 0.45 <= z <= 0.5, {x, y,
     z} ];

I am trying to gave a shot at the solution with following parameters

ba[x_, y_, z_] = 1;
kappa = 4;
sigma = 2;
temp = 0.5;

I tried the following method and got the error:

Needs["NDSolve`FEM`"];

sol = NDSolveValue[{eqns[[1]] == 
     NeumannValue[0, z == 0.45 || z == 0.5], 
    eqns[[2]] == NeumannValue[0, z == 0.45 || z == 0.5], 
    eqns[[3]] == NeumannValue[0, z == 0.45 || z == 0.5], 
    eqns[[4]] == NeumannValue[0, z == 0.45 || z == 0.5] , 
    u3[0, x, y, z] == 0.5*x*ba[ x, y, z], 
    u4[0, x, y, z] == 0.5*x*ba[x, y, z]}, {u1, u2, u3, u4}, {t, 0, 
    10}, {x, y, z} \[Element] region3d, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {MaxCellMeasure -> 0.005}}];

NDSolveValue::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable. >>

Since I am solving such an equation for the first time, I am not particularly sure that I have used the Neumann boundary conditions correcly.

Answer 1

I propose to reproduce the 2D model from the article "Magnetic Flux Penetration in a Mesoscopic Superconductor with a Slit by Isaias G. de Oliveira". We use equations (1), (2) as follows

A = {u3[t, x, y], u4[t, x, y]};
psi = u1[t, x, y] + I u2[t, x, y];
sigma D[A, t] + 
  1/2/k I ((u1[t, x, y] - 
        I u2[t, x, y]) Grad[(u1[t, x, y] + I u2[t, x, y]), {x, 
        y}] - (u1[t, x, y] + 
        I u2[t, x, y]) Grad[(u1[t, x, y] - I u2[t, x, y]), {x, 
        y}]) + (u1[t, x, y]^2 + u2[t, x, y]^2) A - 
  Laplacian[A, {x, y}] // Simplify

D[psi, t] + I/k A.Grad[psi, {x, y}] - 1/k^2 Laplacian[psi, {x, y}] + 
  A.A psi -(1-T) psi (1 - (u1[t, x, y]^2 + u2[t, x, y]^2)) // Simplify

ComplexExpand[%] 

Putting all the equations together we find

    eqns={-u1[t, x, y] + T*u1[t, x, y] + u1[t, x, y]^3 - T*u1[t, x, y]^3 + u1[t, x, y]*u2[t, x, y]^2 - T*u1[t, x, y]*u2[t, x, y]^2 + u1[t, x, y]*u3[t, x, y]^2 + 
  u1[t, x, y]*u4[t, x, y]^2 - (u4[t, x, y]*Derivative[0, 0, 1][u2][t, x, y])/k - Derivative[0, 0, 2][u1][t, x, y]/k^2 - 
  (u3[t, x, y]*Derivative[0, 1, 0][u2][t, x, y])/k - Derivative[0, 2, 0][u1][t, x, y]/k^2 + Derivative[1, 0, 0][u1][t, x, y], 
 -u2[t, x, y] + T*u2[t, x, y] + u1[t, x, y]^2*u2[t, x, y] - T*u1[t, x, y]^2*u2[t, x, y] + u2[t, x, y]^3 - T*u2[t, x, y]^3 + u2[t, x, y]*u3[t, x, y]^2 + 
  u2[t, x, y]*u4[t, x, y]^2 + (u4[t, x, y]*Derivative[0, 0, 1][u1][t, x, y])/k - Derivative[0, 0, 2][u2][t, x, y]/k^2 + 
  (u3[t, x, y]*Derivative[0, 1, 0][u1][t, x, y])/k - Derivative[0, 2, 0][u2][t, x, y]/k^2 + Derivative[1, 0, 0][u2][t, x, y], 
 u1[t, x, y]^2*u3[t, x, y] + u2[t, x, y]^2*u3[t, x, y] - Derivative[0, 0, 2][u3][t, x, y] + (u2[t, x, y]*Derivative[0, 1, 0][u1][t, x, y])/4 - 
  (u1[t, x, y]*Derivative[0, 1, 0][u2][t, x, y])/4 - Derivative[0, 2, 0][u3][t, x, y] + 2*Derivative[1, 0, 0][u3][t, x, y], 
 u1[t, x, y]^2*u4[t, x, y] + u2[t, x, y]^2*u4[t, x, y] + (u2[t, x, y]*Derivative[0, 0, 1][u1][t, x, y])/4 - (u1[t, x, y]*Derivative[0, 0, 1][u2][t, x, y])/4 - 
  Derivative[0, 0, 2][u4][t, x, y] - Derivative[0, 2, 0][u4][t, x, y] + 2*Derivative[1, 0, 0][u4][t, x, y]};

We use the data from the article and find a solution

L = 1; region2d = 
 ImplicitRegion[-L <= x <= L && -L <= y <= 
    L && ! (0 < x <= L && -0.1 < y < 0.1), {x, y}];
ba[x_, y_] := 1;
k = 4;
sigma = 2; T = .5; psi0 = 1/Sqrt[2 ];



ic = {u1[0, x, y] == psi0, u2[0, x, y] == psi0, u3[0, x, y] == 0, 
   u4[0, x, y] == 0};
bc = DirichletCondition[{u1[t, x, y] == psi0, u2[t, x, y] == psi0, 
    u3[t, x, y] == -0.5*y*ba[x, y] (1 - Exp[-5 t]), 
    u4[t, x, y] == 0.5*x*ba[x, y] (1 - Exp[-5 t])}, True];

t0 = 2; sol = 
 NDSolveValue[{eqns == {0, 0, 0, 0}, bc, ic}, {u1, u2, u3, u4}, {t, 0,
    t0}, {x, y} \[Element] region2d];

Table[DensityPlot[sol[[i]][t0, x, y], {x, y} \[Element] region2d, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  PlotPoints -> 40, PlotRange -> All], {i, 4}]

I'm not sure that quantum vorticity can be obtained in this model. The author of the article missed the details. It is known that he used the COMSOL, but how the vortices arise there is not clear.