3
$\begingroup$

I am solving the motion equation of an isolated vortex oscillation in a superconductor. I am assuming that the driving force is generated by an RF current oscillating at 1.3 GHz. The main point of this calculation is to see what's the effect of flux-flow instability under RF drive. The code I wrote is reported below.

f = 1.3 10^9; (*Hz*)
\[Omega] = 2. Pi f;
Trf = 1./f;

Tc = 9.25; (*K*)
T = 1.5; (*K*)

\[Mu]0 = 4. Pi 10^-7; (*H/m*)
\[Phi]0 = 2.07 10^-15; (*Wb*)
n = 5.56 10^28; (*m^-3*)
e = 1.6 10^-19; (*C*)
m = 9.1 10^-31; (*kg*)

\[Tau] = 0.5 10^-9; (*s*)
vf = 1.37 10^6; (*m/s*)
Bc2 = 410. 10^-3; (*T*)
l = 100 10^-9; (*m*)
\[Lambda] = 49.5 10^-9; (*m*)
\[Xi] = 28 10^-9; (*m*)
\[Kappa] = \[Lambda]/\[Xi];

\[Gamma] = (\[Phi]0)/(\[Mu]0 \[Lambda] );

v0 = Sqrt[(l vf Sqrt[14 Zeta[3] (1 - T/Tc)])/(3 Pi \[Tau])];

\[Rho] = (m vf)/(n e^2 l);
\[Eta] = (\[Phi]0 Bc2)/\[Rho];

g = Log[\[Kappa]]+0.5+Exp[-0.4-0.8 Log[\[Kappa]]-0.1(Log[\[Kappa]])^2]; 
\[Epsilon] = (\[Phi]0^2 g)/(4. Pi \[Mu]0 \[Lambda]^2);

Tmax = 3. Trf;
Zmax = 1. 10^-6; (*m*)

B = 100. 10^-3; (*T*)

sol = NDSolve[{
     (v0^2 \[Eta])/(v0^2 + D[x[t, z], t]^2)D[x[t, z],t] == \[Epsilon] D[x[t, z], z, z] + \[Gamma] B Cos[\[Omega] t] Exp[-(z/\[Lambda])],
     x[0., z] == 0.,
     x[t, Zmax] == 0.,
     (D[x[t, z], z] /. z -> 0) == 0.
     }, x, {t, 0., Tmax}, {z, 0., Zmax},
    AccuracyGoal -> 13,
    PrecisionGoal -> 2,
    MaxSteps -> Infinity];
u[t_, z_] = Evaluate[x[t, z] /. sol][[1]];
du[t_, z_] = D[u[t, z],t];

phaseSpace = 
   ParallelTable[{Re[u[t, 0.]] 10^6, Re[du[t, 0.]] 10^-3}, {t, Trf, 
     Tmax, Tmax/500.}];
ListPlot[phaseSpace, Joined -> True, PlotRange -> All]

I am able to get a almost correct result out of this calculation only if I fix the solution domain to 1 micron (Zmax = 1. 10^-6). This is however not correct because the vortex should be able to oscillate more than that inside the bulk, therefore I am literally chopping my solution using a too small domain.

If I fix Zmax to a larger number (e.g. 5 um) I need to fix MaxStepSize to a small number in order to avoid errors and the calculation would take more than 2 days on my PC (I actually stopped it after 2 days because the PC froze up!). Does anybody know how to solve this problem?

Also, I was wondering if there is any way to solve the equation when the vortex velocity is higher than v0 (the velocity onset for flux-flow instability). If I fix B to let's say 200 mT, then the vortex velocity would be larger than v0 and the calculation does not return a correct solution, just a series of spikes. Is there any Method I can use to solve this problem?

Please let me know! Thank you in advance.

Mattia

Improve this question
$\endgroup$
2
  • $\begingroup$ I posted an answer, but it is probably a numerical artifact. I noticed that all of your BCs and ICs are zero. Is that correct? $\endgroup$
    – Tim Laska
    Commented Jun 13, 2019 at 10:55
  • $\begingroup$ Yes it is correct. In the bulk the vortex do not move and it is perpendicular to the surface. However, the initial condition might be set differently I believe. The best condition would be to put the IC for a vortex with zero velocity, but I do not know the analytical form of it, that's why I put zero as IC (static vortex). I figured that if the solution is calculated for a large enough number of oscillations, it converges to a steady state solution anyway. $\endgroup$
    – Mattia
    Commented Jun 17, 2019 at 15:46

2 Answers 2

Reset to default
5
$\begingroup$

When working with small spatial and temporal scales, it is often conducive to scale the variables by characteristic lengths and times (or reciprocal frequencies). The natural choices for length is $\lambda$ and reciprocal time is $\omega$. If we create dimensionless $x$, $z$, and $t$, and do the appropriate substitutions, we should be able to derive a dimensionless differential equation that depends on several dimensionless groups.

My non-dimensionalized version of your dimensioned equation is

$$\frac{{\frac{{\partial x}}{{\partial t}}}}{{{\nu ^2}{{\left( {\frac{{\partial x}}{{\partial t}}} \right)}^2} + 1}} = \delta \left( {\rho \frac{{{\partial ^2}x}}{{\partial {z^2}}} + Cos(t)Exp( - z)} \right)$$

If we use your initial parameters:

f = 1.3 10^9;(*Hz*)ω = 2. Pi f;
Trf = 1./f;

Tc = 9.25;(*K*)
T = 1.5;(*K*)
μ0 = 4. Pi 10^-7;(*H/m*)
ϕ0 = 2.07 10^-15;(*Wb*)
n = 5.56 10^28;(*m^-3*)
e = 1.6 10^-19;(*C*)
m = 9.1 10^-31;(*kg*)
τ = 0.5 10^-9;(*s*)
vf = 1.37 10^6;(*m/s*)
Bc2 = 410. 10^-3;(*T*)
l = 100 10^-9;(*m*)
λ = 49.5 10^-9;(*m*)
ξ = 28 10^-9;(*m*)
κ = λ/ξ;

γ = (ϕ0)/(μ0 λ);

v0 = Sqrt[(l vf Sqrt[14 Zeta[3] (1 - T/Tc)])/(3 Pi τ)];

ρ = (m vf)/(n e^2 l);
η = (ϕ0 Bc2)/ρ;

g = Log[κ] + 0.5 + 
   Exp[-0.4 - 0.8 Log[κ] - 0.1 (Log[κ])^2];
ϵ = (ϕ0^2 g)/(4. Pi μ0 λ^2);

Tmax = 3. Trf;
Zmax = 1. 10^-6;(*m*)
B = 100. 10^-3(*T*)

Here is my modified non-dimensionalized version of your code.

Zmax = 100;
Tmax = 20;
vc = λ ω;
nu = vc/v0;
rho = ϵ/(λ B γ);
delta = (B γ)/(η vc);
sol = NDSolve[{D[x[t, z], t]/(1 + nu^2 D[x[t, z], t]^2) == 
     delta ( rho D[x[t, z], z, z] + Cos[ t] Exp[-z]), x[0., z] == 0., 
    x[t, Zmax] == 0., (D[x[t, z], z] /. z -> 0) == 0.}, 
   x, {t, 0., Tmax}, {z, 0., Zmax}];
u[t_, z_] = Evaluate[x[t, z] /. sol][[1]];
du[t_, z_] = D[u[t, z], t];
phaseSpace = 
  Table[{Re[u[t, 0.]] , Re[du[t, 0.]] }, {t, Trf ω, Tmax, 
    Tmax/500.}];
ListPlot[phaseSpace, Joined -> True, PlotRange -> All]

Phase plot

It took about 2.5 seconds to run on my machine. I won't have time to double check my work for a couple days, but I hope this gets you started.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I received two messages: NDSolve::mxst: Maximum number of 10000 steps reached at the point t == 2.888363945453023*^-45.` NDSolve::eerr: Warning: scaled local spatial error estimate of 576279.2946597305 at t = 2.888363945453023*^-45 in the direction of independent variable z is much greater than the prescribed error tolerance. $\endgroup$
    – Alex Trounev
    Commented Jun 13, 2019 at 4:04
  • $\begingroup$ You are correct. I did not hear the beeps last night, but I hear them now and I see the same message. More investigation is needed. $\endgroup$
    – Tim Laska
    Commented Jun 13, 2019 at 10:10
  • $\begingroup$ @AlexTrounev I noticed that all the BCs and ICs are zero, which seems like it would be difficult/impossible to evolve. I'll ask the OP for additional guidance. $\endgroup$
    – Tim Laska
    Commented Jun 13, 2019 at 10:51
  • $\begingroup$ Thank you for the suggestion to make everything dimensionless, I believe it helped a lot on speeding up the calculation. I get the same error Alex Trounev get, but the code still spits out the solution. $\endgroup$
    – Mattia
    Commented Jun 17, 2019 at 16:07
2
$\begingroup$

A small code change allows you to solve both problems. Figure 1 shows solutions for increasing B and Zmax.

f = 
 1.3 10^9;(*Hz*)\[Omega] = 2. Pi f;
Trf = 1./f;

Tc = 9.25;(*K*)T = 1.5;(*K*)\[Mu]0 = 4. Pi 10^-7;(*H/m*)\[Phi]0 = 
 2.07 10^-15;(*Wb*)n = 5.56 10^28;(*m^-3*)e = 1.6 10^-19;(*C*)m = 
 9.1 10^-31;(*kg*)\[Tau] = 0.5 10^-9;(*s*)vf = 1.37 10^6;(*m/s*)Bc2 = 
 410. 10^-3;(*T*)l = 100 10^-9;(*m*)\[Lambda] = 
 49.5 10^-9;(*m*)\[Xi] = 28 10^-9;(*m*)\[Kappa] = \[Lambda]/\[Xi];

\[Gamma] = (\[Phi]0)/(\[Mu]0 \[Lambda]);

v0 = Sqrt[(l vf Sqrt[14 Zeta[3] (1 - T/Tc)])/(3 Pi \[Tau])];

\[Rho] = (m vf)/(n e^2 l);
\[Eta] = (\[Phi]0 Bc2)/\[Rho];

g = Log[\[Kappa]] + 0.5 + 
   Exp[-0.4 - 0.8 Log[\[Kappa]] - 0.1 (Log[\[Kappa]])^2];
\[Epsilon] = (\[Phi]0^2 g)/(4. Pi \[Mu]0 \[Lambda]^2);

Tmax = 3. Trf;
Zmax = 5. 10^-6;(*m*)B = 300. 10^-3;(*T*)sol = 
 NDSolveValue[{ \[Eta]/(1 + D[x[t, z], t]^2/v0^2) D[x[t, z], 
       t] == \[Epsilon] D[x[t, z], z, 
        z] + \[Gamma] B Cos[\[Omega] t] Exp[-(z/\[Lambda])], 
    x[0., z] == 0., 
    x[t, Zmax] == 0., (D[x[t, z], z] /. z -> 0) == 
     0.}, {Re[x[t, 0]]*10^6, Re[D[x[t, 0], t]]*10^-3}, {t, 0., 
    Tmax}, {z, 0., Zmax}, AccuracyGoal -> 13, PrecisionGoal -> 2, 
   MaxSteps -> Infinity] // Quiet;


phaseSpace = ParallelTable[sol, {t, Trf, Tmax, Tmax/500.}];
ListPlot[phaseSpace, Joined -> True, PlotRange -> All, 
 PlotLabel -> Grid[{{"Zmax = ", Zmax}, {"B = ", B}}]]

Figure 1

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you Alex. There is something strange in the solutions you get, as a function of B both speed (y-axis) and displacement (x-axis) should increase as a function of the rf field amplitude B. I noticed as well that if Zmax is increased (keeping the same B) the code tends to return a solution that qualitatively is ok (an ellipse) but it is not accurate and generally underestimated. To partially solve this I used the "SingleStep" method specifying MaxStepSize and StartingStepSize. The solution in this way is more accurate when Zmax is increased. $\endgroup$
    – Mattia
    Commented Jun 17, 2019 at 16:11
  • $\begingroup$ I meant "FixedStep", sorry $\endgroup$
    – Mattia
    Commented Jun 17, 2019 at 16:17
  • $\begingroup$ @Mattia OK! You probably understand what the solution should be. I just indicated how to choose a branch of solution using ` NDSolveValue` instead of ` NDSolve`. $\endgroup$
    – Alex Trounev
    Commented Jun 17, 2019 at 16:36

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.