Hello,
I would like to calculate the current distribution in a set of (superconducting) striplines (they are parallel and “infinite”).
Usually, the total current I_bias circulating in a (single) stripline is an input data/constraint for this problem (while the distribution in the stripline is far from being uniform). How would it be possible to implement such constraint ( iint from{Supercond} u (x,y)dS = I_bias ) with the basic installation of FreeFEM ?
or do I need to install additional software (IPOPT ?), and possibly recompile the code ?
Thank you for any advice.
Denis
I’m not sure I understand you correctly, but can’t you discretize your constraint using a varf in a vector c, and then solve this using a Lagrangian, e.g., matrix B = [[A,c],[c',0]];? Here is such an example where we impose a mean of 0. You can strip the PETSc stuff and use plain FreeFEM matrices.
Thank you very much for your suggestion. I have made quite a number a attempts to understand how to manipulate the stiffness matrix, the Lagrange function, and even the boundary conditions at the edge of the space domain. I may have to return to the basics and provide more details:
My objective is to calculate the current density in a system of superconducting wires. Input data are the geometry and the total currents circulating in the different wires. As in the general case the system does not have any symmetry (except translational invariance along direction z), I chose not to use cylindrical coordinates.
The script below simulates a single cylindrical infinite wire made with superconducting material and displays the following graphs: mesh, vector-potential component Az, and current density Jz in the wire.
For this particular example, the solution is well-known, Az ~ log® in the “air” (= outside the wire). Thus I can use the boundary condition Az=mu_0I/4pilog(x^2 +y^2) on the wall limiting the domain. Line 22 of the script can be modified to
-
change the space between rectangular section and circular sections of the box
-
change the script usng “PROBLEM/SOLVE” or using the stiffness matrix A : we obtain very similar results.
Meissner-supra_wire-3.edp (4.7 KB)
=====BEGINNING OF REPORT (Meissner-supra_wire-3.edp) ===============
– FreeFem++ v 3.380001 (date Sat Feb 6 20:00:25 UTC 2016)
Load: lg_fem lg_mesh lg_mesh3 eigenvalue
1 : // Vector potential created by an infinite superconducting wire
2 : // Solved using either the “problem” keyword, or manipulating the stiffness matrix.
3 : // Boundary conditions at large distance being non-zero, a current must circulate
4 : // Where Laplace(Az) is not 0, i.e. in the wire
5 : // “boxoice” parameter controls whether we use a space with rectangular or cylindrical section.
6 : ///////////////Define geometry///////////////////////////////
7 : // linear dimensions (in µm)
8 : // Penetration depth (squared)
9 : real tgv=1E30, epsi=1E-6, xi,yi, Lond=0.2, l2= Lond^2; // Lambda_London =0.2 um
10 : real mu0s2p = 2e-7, mu0=mu0s2p2pi, I1=1, A0=mu0s2pI1; // for current I1=1 Amp.
11 : // Rectangular Section: Choice 1 = Results are marginally acceptable
12 : real far=20,high=20;
13 : border bottom(t=-far,far) { x=t; y=-high; label=1; }
14 : border rside(t=-high,high) { x=far; y=t; label=1; }
15 : border lside(t=high,-high) { x=-far; y=t; label=1; }
16 : border top(t=far,-far) { x=t; y=high; label=1; }
17 :
18 : // Cylindrical section: Choice 2 = Results are just nonsense…
19 : real RB=23; // RB is the radius of the box
20 : border RoundBox(t=0,2pi) { x=RBcos(t); y=RBsin(t); label=2; }
21 : Meissner-supra_wire_5bi4.edp
22 : int boxoice= 2, prost=1; // boxoice: choice of outer box shape; prost: 0 for “problem”; 1 for “stiffness”
23 :
24 :
25 : // Wire
26 : real R=1; // R is the radius of the wire
27 : border wire(t=0,2pi) { x=Rcos(t); y=Rsin(t); label=3; }
28 : border bulk(t=0,2pi) { x=(R-2Lond)cos(t); y=(R-2Lond)sin(t); label=4; }
29 : /////////////////////////////////////////////////////////////
30 : //////////// Mesh and London equation ///////////////////////
31 : int NBM=50, NBW=200;
32 : mesh Th;
33 : if ( boxoice==1 Meissner-supra_wire_5bi4.edpMeissner-supra_wire_5bi4.edp) Th=buildmesh(bottom(NBM)+rside(NBM)+top(NBM)+lside(NBM)+wire(NBW)+bulk(NBW));
34 : if ( boxoice==2 ) Th=buildmesh(RoundBox(3NBM)+wire(NBW)+bulk(NBW));
35 : // regions are numbered 0 , 1 and 2
36 : int core=Th(0,0).region, super=Th(R-Lond,0).region, air=Th(R+Lond,0).region;
37 : cout << “Superc. regions: " << core <<” and "<< super << ", and air: "<< air<< endl;
38 :
39 : fespace Vh(Th,P2); // With P1, the result would be very noisy…
40 : Vh Az, v, u0, ubc;
41 : int nbvertices=Th.nv, NDF=Vh.ndof;
42 :
43 : // The following lines give analytically exact BC for single cylindrical wire
44 : if ( boxoice==1 ) {
45 : u0 =(abs(abs(x)-far)<1E-2 || abs(abs(y)-high)<1E-2) ? A0/2log(x^2 +y^2):0; // for rectangular box
46 : };
47 : if ( boxoice==2 ) {
48 : u0 =(abs(x^2+y^2-RB^2)<epsi) ? A0log(RB):0; // for disk box
49 : };
50 : ubc=(u0==0) ? 0: tgvu0; // This was to avoid 01E30, but it does not seem necessary …
51 :
52 : fespace W0(Th,P0);
53 : W0 Sup=(region==super || region==core); // Size(Sup)=nb of TRIANGLES !!!
54 : W0 Cor=(region==core);
55 : cout << “Nb Vert.=”<< nbvertices<<", NDoF="<< NDF<< “, size(Sup)=” << Sup.n << endl;
56 : func fL=-Sup/l2; // For London equation in the superconductor
57 : ///////////////////////////////////////////////////
58 : //////////// Solve problem ///////////////////////
59 : problem poissonMagneto(Az,v) =
60 : int2d(Th)( dx(Az)dx(v) + dy(Az)dy(v) )
61 : + int2d(Th) ( -vAzfL ) // the source term: Mu0J=-Az/Lond^2 in the superconductor
62 : + on(boxoice, Az=u0); // boundary condition far away: Az=mu0s2pI1log®
63 :
64 : varf va(Az,v) = // definition of the problem
65 : int2d(Th)( dx(Az)dx(v) + dy(Az)dy(v))
66 : + int2d(Th) ( -vAzfL ) // the source term: Mu0J=-Az/l2 in the superconductor
67 : + on(boxoice,Az=u0);
68 : matrix A=va(Vh,Vh);
69 : real[int] F=ubc[];
70 :
71 : if (prost==0) {poissonMagneto;}
72 : else {
73 : set(A,solver=UMFPACK);
74 : Az[]=A^-1F; // solve the linear system
75 : }; // prost==1 seems to be slightly faster; otherwise, both methods give the same results.
76 :
77 : ////////////////////////////////////////////////////
78 : //////////// Display results ///////////////////////
79 :
80 : plot(Th,ps=“wire.eps”,bw=1, wait=true, cmm=“Mesh”);
81 :
82 : Vh Bx = -dy(Az)(x,y), By = dx(Az)(x,y);
83 : Vh Jz = (dx(By)-dy(Bx))/mu0;
84 : string Title="Contour plot of Az, ";
85 : if (prost==0) Title=Title+“Solve”; else Title=Title+“Stiffness”;
86 : plot(Az, ps=“vector_potential_wire.eps”, value=true, wait=true, cmm=Title);
87 : Title=“Contour plot of Jz, “;
88 : if (prost==0) Title=Title+“Solve”; else Title=Title+“Stiffness”;
89 : plot(Jz, ps=“Current_density_wire.eps”, value=true, wait=true, cmm=Title);
90 : // Larger valuesMeissner-supra_wire_5bi4.edp of the current are indeed obtained for the surface of the wire,
91 : cout << “Core current =” << int2d(Th)(CorJz) << endl; // 9% of the current in the core
92 : // No negative values in the core, as expected
93 : // Check total current: integrates the current density in the wire =I1
94 : cout << “Total current =” << int2d(Th)(Jz) << " and Wire current =" << int2d(Th)(SupJz) << endl;
95 : // The integrated current density is about 10% too small
96 :
97 : // Using Gauss theorem, calculate By(R,0)=mu0s2pI1/R, i.e. mu0s2pI1=RBy(R,0)
98 : cout << "Apply Gauss theorem at r=2R: I1=”<< 2RBy(2R,0)/mu0s2p << endl;
99 : // The value of the field is also too small (about -7%). Does not seem to depend on the mesh… ???
100 : cout << "Apply Gauss theorem at r=20R: I1=”<< 20RBy(20*R,0)/mu0s2p << endl;
101 : // Here also, By is too small (-9%).
102 :
103 : sizestack + 1024 =7872 ( 6848 )– mesh: Nb of Triangles = 42176, Nb of Vertices 21164
Superc. regions: 0 and 1, and air: 2
Nb Vert.=21164, NDoF=84503, size(Sup)=42176
Core current =0.0295082
Total current =83.7838 and Wire current =0.29396
Apply Gauss theorem at r=2R: I1=0.302521
Apply Gauss theorem at r=20R: I1=0.302628
times: compile 0.018529s, execution 5.9965s, mpirank:0
CodeAlloc : nb ptr 3284, size :368096 mpirank: 0
Ok: Normal End
===========END OF REPORT==================
This script gives unreliable results for circular section of the box limiting the space. When the space is defined with a rectangular section, the results are much better but not fully reliable:
- Larger values of the current are indeed obtained for the surface of the wire,
- No negative values in the core, as expected
- The integrated current density is about 10% too small
- the magnetic field “By” roughly follows a 1/r law.
- The value of “By” is also too small (about -7%). Does not seem to depend on the mesh… ???
How is it possible to improve the results ?
The approaches presented above allow in principle to compute the current circulating in the wire. However, with a more complex system of conductors, the boundary conditions become difficult to evaluate beforehand. Is it possible to impose the value of the current I1 circulating in the wire (int2d(Th2, Jz) =I1) and calculate Az without giving boundary conditions on borders labeled “1” (or “2” in the case of the cylinbdrical box)?
I guess that it is necessary to manipulate the stiffness matrix A and the vector F associated to the linear form to implement the condition on the total currents circulating in the wires.
If I understand correctly, the equation A.u=F is not 100% correct because of edge effects: the points belonging to the border are missing neighbours to evaluate the 2nd derivative along the axis which is perpendicular to the border. This is presumably the fundamental reason why you must change the coefficients of the stiffness matrix for lines corresponding to points on the border, for a set of proper boundary conditions.
The second code (Meissner-supra_wire-5bi4.edp) attempts to implement the condition for the total current in a wire, using the algorithm of the Lagrange multiplier. This can be summarized as follows:
-find the solution of \Delta u = fLu given that \int u dS = mu0\London^2I, with fL=1/London^2 in the superconducting wire, and 0 elsewhere. For c = vector of triangle areas, it is equivalent to solve the matrix equations Au=F and c’u= II=mu0\Lambda^2*I, which result from the minimization of the Lagrage \int (U^2/2 -fL) dS - Lambda (\int u dS -I), where lambda is the Lagrange multiplier.
Therefore, I calculated the solution of AA * u= [[A, c];[c’,0]]U = [ZERO,II], where ZERO is the presumably the vector of NDoF elements, each of them equal to 0, and U = [u,Lambda] : Au+cLambda=ZERO; c’*u = II.
Because not giving boundary conditions is equivalent to giving Neumann boundary conditions, I still must include realistic BC, i.e. I must change the vector ZERO elements corresponding to the border. In a first try, I just want to give the Dirichlet boundary conditions corresponding to the expected analytical result, i.e. the same as in Meissner-supra_wire-3.edp.
The second code has 2 toggle parameters, boxoice and solvoice, the first one allows to choose between a rectangluar section box or a disk section box; the second parameter is used to choose either to solve the matrix exactly as in Meissner-supra_wire-3.edp, or to solve by minimization of the Lagrange function.
Meissner-supra_wire_5bi4.edp (6.7 KB)
========== BEGINNING OF REPORT (Meissner-supra_wire_5bi4.edp) =================
– FreeFem++ v 3.380001 (date Sat Feb 6 20:00:25 UTC 2016)
Load: lg_fem lg_mesh lg_mesh3 eigenvalue
1 : // Vector potential created by an infinite superconducting wire (circular cross-section)
2 : // Attempt to implement the total current condition with the “augmented matrix” approach
3 : // “boxoice” parameter controls whether we use a rectangular space or a disk.
4 :
5 : ///////////////Define geometry///////////////////////////////
6 : // linear dimensions (in µm)
7 : // Penetration depth (squared)
8 : real tgv=1E30, epsi=1E-10, xi,yi,Lond=0.2, l2= Lond^2; // Lambda_London =0.2 um
9 : real mu0s2p = 2e-7, mu0=mu0s2p2pi, I1=1, A0=mu0s2pI1; // for current I1=1 Amp.
10 :
11 : // Rectangular Box: Choice 1
12 : real far=20,high=20;
13 : border bottom(t=-far,far) { x=t; y=-high; label=1; }
14 : border rside(t=-high,high) { x=far; y=t; label=1; }
15 : border lside(t=high,-high) { x=-far; y=t; label=1; }
16 : border top(t=far,-far) { x=t; y=high; label=1; }
17 :
18 : // Disk box: Choice 2
19 : real RB=23; // RB is the radius of the box
20 : border RoundBox(t=0,2pi) { x=RBcos(t); y=RBsin(t); label=2; }
21 :
22 : int boxoice= 2, solvoice=0; // boxoice: choice of outer box shape;
23 : // solvoice : choice of matrix. 0:modified stiffness matrix; 1:original stiffness matrix
24 :
25 : // Superconducting wire
26 : real R=1; // R is the radius of the superconducting wire
27 : border wire(t=0,2pi) { x=Rcos(t); y=Rsin(t); label=3; }
28 : border bulk(t=0,2pi) { x=(R-2Lond)cos(t); y=(R-2Lond)sin(t); label=4; }
29 : /////////////////////////////////////////////////////////////
30 : //////////// Mesh and London equation ///////////////////////
31 : int NBM=50, NBW=250;
32 : mesh Th;
33 : if ( boxoice==1 ) Th=buildmesh(bottom(NBM)+rside(NBM)+top(NBM)+lside(NBM)+wire(NBW)+bulk(NBW));
34 : if ( boxoice==2 ) Th=buildmesh(RoundBox(3NBM)+wire(NBW)+bulk(NBW));
35 : // regions are numbered 0 , 1 and 2
36 : int core=Th(0,0).region, super=Th(R-Lond,0).region, air=Th(R+Lond,0).region;
37 : cout << “Superc. regions: " << core <<” and "<< super << ", and air: "<< air<< endl;
38 :
39 : fespace Vh(Th,P2); // P1 gives too much noise.
40 : Vh Az, v, u0, ubc, Sup=(region==super || region==core); // Size(Sup)=Vh.ndof (nb of TRIANGLES with P0 and P1);
41 : int NDF=Vh.ndof;
42 : fespace W0(Th,P0);
43 : W0 Cor=(region==core);
44 : varf varea(unused,chiK) =int2d(Th)(chiK); // the space function constant / triangle
45 : real[int] c=varea(0,Vh);
46 : // I want to zero the elements of c not associated to a vertex (or triangle) in the superconductor
47 : // c=c.Sup[]; // Wrong: results do not look like results of Meissner-supra_wire-3.edp
48 : // I don’t understand why there is such a difference in integrating the current density only on the wire,
49 : // compared to integrating on the full space, i.e. on the wire + regions where the current density should be 0!!
50 : cout << “Size(Sup)=” << Sup.n << ", min© = " << c.min << ", max© = " << c.max << ", sum© = " << c.sum <<endl;
51 :
52 : func fL=-Sup/l2; // For London equation in the superconductor
53 :
54 : ///////////////////////////////////////////////////
55 : //////////// Solve problem ///////////////////////
56 :
57 : real lm;
58 : // Is “Not giving boundary conditions” equivalent to “Neumann” ?
59 : // It is OK for the case with cylindical symmetry, but not in the general case
60 : // In the rectangular space, do we need to give something more realistic ?
61 : // Or alternatively, how can we suppress the default Neumann condition ?
62 : // It should be enough to specify the source term with “I1”…
63 : // The following lines give analytically exact BC for single cylindrical wire
64 : if ( boxoice==1 ) {
65 : u0 =(abs(abs(x)-far)<1E-2 || abs(abs(y)-high)<1E-2) ? A0/2log(x^2 +y^2):0; // for rectangular box
66 : ubc=(abs(abs(x)-far)<1E-2 || abs(abs(y)-high)<1E-2) ? tgvA0/2log(x^2 +y^2):0; // for BC
67 : };
68 : if ( boxoice==2 ) {
69 : u0 =(abs(x^2+y^2-RB^2)<epsi) ? A0log(RB):0; // for disk box
70 : ubc=(abs(x^2+y^2-RB^2)<epsi) ? tgvA0log(RB):0; // for BC
71 : };
72 : real[int] F = ubc[], FI=[F,-mu0l2I1], solu(NDF+1); // The source terms are all contained in FI +F (non0)
73 : cout << “NDoF=”<< NDF<<", size© = " << c.n << ", size(u0) = " << u0.n << endl;
74 :
75 : varf va(Az,vh) = // definition of the problem
76 : int2d(Th)( dx(Az)dx(vh) + dy(Az)dy(vh))
77 : + int2d(Th) ( -vhAzfL ) // the source term: Mu0J=-Az/l2 in the superconductor
78 : + on(boxoice,Az=u0) // applies the BC to the proper border
79 : ;
80 :
81 : matrix A=va(Vh,Vh), AA = [ [ A , c] , // As for Lagrange multiplier
82 : [ c’, 0] ] ; // “Total current” constraint: area’Az =Imu0s2p
83 :
84 : if (solvoice==0) {
85 : set(AA,solver=UMFPACK);
86 : solu=AA^-1FI; // solve the linear system; BC are not correct at least for rectangular box
87 : // Instead of having iso-Az looking like circles, they gradually look like rectangles
88 : // when they approach the boundary of the domain, as if the whole boundary was set to the same value
89 : [Az[], lm] = solu;
90 : }
91 : else {
92 : set(A,solver=UMFPACK);
93 : Az[]=A^-1F; // solve the linear system: ONLY RECTANGULAR BOX gives results which are
94 : // reasonable = identical to what is found in Meissner-supra_wire-3.edp
95 : };
96 : ////////////////////////////////////////////////////
97 : //////////// Display results ///////////////////////
98 : // Retrieve the content of A and AA just for display
99 : int[int] I(1),J(1);
100 : real[int] C(1);
101 : [I,J,C]=A;
102 : cout <<“Size(A)=” << A.n << “, size(F)=” << F.n << “, Max(A)=”<< C.max<< “, Min(A)=”<< C.min<< endl;
103 : cout <<“Size(AA)=” << AA.n << “, Max(AA)=”<< C.max<< “, Min(AA)=”<< C.min<< endl;
104 : real MaxF=0, MinF=100, Fi;
105 : for (int i=0; i<F.n; i++) {
106 : if (F(i)>0) {
107 : Fi=F(i)/tgv;
108 : MaxF=MaxF>Fi? MaxF:Fi;
109 : MinF=MinF<Fi? MinF:Fi;
110 : };
111 : };
112 : cout <<"Max(F) = " << MaxF << “, meen(F)=”<< MinF<< “, tgv= “<< tgv<< endl;
113 : plot(Th,ps=“wire.eps”,bw=1, wait=true,cmm=“Mesh Box:”+NBM+” Wire:”+NBW);
114 :
115 : Vh Bx = -dy(Az)(x,y), By = dx(Az)(x,y);
116 : Vh Jz = (dx(By)-dy(Bx))/mu0;
117 : plot(Az, value=true, wait=true, cmm=“Contour plot of Az”);
118 :
119 : plot(Jz, value=true, wait=true, cmm=“Contour plot of Jz”);
120 : // Larger values of the current are indeed obtained for the surface of the wire,
121 : cout << “Core current =” << int2d(Th)(CorJz) << endl; //
122 : // this is ~10% of the current in the surface, same ratio as in “Meissner-supra_wire-3.edp” 
123 : // No negative values in the core, as expected
124 :
125 : cout << "Lagrange factor : " << lm << “, size(Jz) = " << Jz.n << endl;
126 : // Check total current: integrates the current density in the wire =I1 ?
127 : cout << “Total current =” << int2d(Th)(Jz) << " and Wire current =” << int2d(Th)(SupJz) << endl;
128 : // => I do not get what I expect, the current in the wire not equal to I1
129 : // The integrated current density is by far too small. What is in the c matrix ???
130 :
131 : // Using Gauss theorem, calculate By(R,0)=mu0s2pI1/R, i.e. mu0s2pI1=RBy(R,0)
132 : cout << "Apply Gauss theorem at r=2R: I1="<< 2RBy(2R,0)/mu0s2p << endl;
133 : cout << "Apply Gauss theorem at r=4R: I1="<< 4RBy(4R,0)/mu0s2p << endl;
134 : cout << "Apply Gauss theorem at r=20R: I1="<< 20RBy(20R,0)/mu0s2p << endl;
135 : // The above results show whether By(r,0) decreases as 1/r or not.
136 : // only in the central region of the domain for both rectangular and disk boxes
137 : // Is this problem related to the Neumann default boundary conditions ?
138 :
139 : sizestack + 1024 =11400 ( 10376 )
– mesh: Nb of Triangles = 53254, Nb of Vertices 26703
Superc. regions: 0 and 1, and air: 2
Size(Sup)=106659, min© = -1.14275e-16, max© = 0.300256, sum© = 1661.42
NDoF=106659, size© = 106659, size(u0) = 106659
– Block Matrix NxM = 2x2 nxm =106660x106660 nb none zero coef. 1429256
Size(A)=106659, size(F)=106659, Max(A)=1e+30, Min(A)=-2.48824
Size(AA)=106660, Max(AA)=1e+30, Min(AA)=-2.48824
Max(F) = 6.27099e-07, meen(F)=6.27099e-07, tgv= 1e+30
Core current =-0.0418992
Lagrange factor : 3.72496e-09, size(Jz) = 106659
Total current =88.6296 and Wire current =-0.412868
Apply Gauss theorem at r=2R: I1=-0.391594
Apply Gauss theorem at r=4R: I1=-0.279846
Apply Gauss theorem at r=20*R: I1=3.2961
times: compile 0.018462s, execution 11.1552s, mpirank:0
CodeAlloc : nb ptr 3522, size :372400 mpirank: 0
Ok: Normal End
====================END OF REPORT===============
- When the space is “rectangular” and solving for the “original matrix” A (boxoice= 1, solvoice=1), then the calculations are essentially identical to the cases (boxoice= 1, prost=0 or 1) in Meissner-supra_wire-3.edp. And the results are very comparable; it means that results are about 7…9% smaller than expected.
- When the space is “rectangular” and solving for the “modified matrix” AA (boxoice= 1, solvoice=0), then the result become strange, as instead of having iso-Az looking like circles, they gradually look like rectangles when they approach the boundary of the domain, as if the whole boundary was set to the same value. I checked the values stored in F, and they look correct. Taking a look at the currents, the current density in the wire turned to NEGATIVE, which is against all expectations… and spikes of large positive currents show up at the surface of the wire. Total currents integrated over the different regions hardly depend on the value of the last element of FI (expected to fix the current in the wire).
- When the space is cylindrical and solving for the “original matrix” A (boxoice= 2, solvoice=1), then the calculations are essentially identical to the cases (boxoice= 2, prost=0 or 1) in Meissner-supra_wire-3.edp. And the results are very comparable; it means that results are far off…
- When the space is cylindrical and solving for the “modified matrix” AA (boxoice= 2, solvoice=0), then the results are as strange as in Meissner-supra_wire-3.edp; it means that results are far off, currents in the wire are NEGATIVE… etc. And the total current is not even comparable to the total current calculated with Meissner-supra_wire-3.edp, cylindrical cases.
I am convinced I missed something in the manipulation of the stiffness matrix, but I am really surprised to see that the shape of the space matters!
Thank you very much for any help.
Denis
I know that this reply is a bit late, but because I am also trying to solve some superconducting geometries I came across your post. You might want to look at the paper “Meissner-London state in superconductors of rectangular cross section in a perpendicular magnetic field”, PHYSICAL REVIEW B 62, 115 (2000). They solve for the response of a long superconducting strip in a magnetic field. The equation they solve is ∇A -A/λ^2, with the penetration length λ = 0 outside the strip. Here, the vector potential A has only a z-component (where z is perpendicular to the 2D plane of the simulation, i.e., z is parallel to the strip). This problem is easy to solve in FreeFem.
The boundary condition they use is A = -B0*x at the top and bottom of the boundary, to give a uniform field far from the strip. But I have also solved this problem for a fixed current I in the strip as follows. We know that if a strip has a current I, then far from the strip it will look like a circular wire and A = constant at a distant value of r. So you can make a circular boundary and set A = constant, and then adjust this constant so that the integrated current in the wire is I.
Because the London equations are linear, I think you can do this for each strip to find the current distribution inside the strip due to I, and the field outside it. I’m not sure how to put the various strips next to each other, though.
Hello,
thank you very much for the interest in my problem. I’ll try to model a superconducting pellet as done in the paper you referenced:
This will probably help me in checking that the code I developed gives comparable results.
I agree that the potential A=constant at a large distance r; more precisely, this constant varies as log(r). What I don’t understand is the difference in taking a circular Boundary with the condition A=constant, or, a rectangular boundary with A~log(x²+y²): the latter case breaks the cylindrical symmetry of the potential, even when the wire is indeed cylindrical!
Referring to the case of several superconducting strips, I don’t know either how to impose currents in each of the conductors separately. If we have 2 conductors and feed +I in the first conductor and -I in the second, the vector potential at long distance is 0 (first order) and trying to apply non-cylindrical boundary is as helpless as in the situation described in the above paragraph (with rectangular boundary). In the best case, it gives identically 0 for the current distribution.
I have not been able to model the influence of a current circulating in one of the wire inducing screening currents in the other superconductors. I am still very much interested in finding a solution, moreover by these times where I cannot access other licensed commercial software.
Thank you very much again.
Best regards
Denis
I thought a bit more about this. Consider the case you mentioned with two strips, with +I and -I. I think you could solve this as follows. First, find the vector potential A1(x, y) due to the current one strip by the method I outlined earlier (surround it by a circular boundary with A = constant). A1 needs to be saved, perhaps by exporting it to disk.
Then, separately, do the same for the second strip, finding A2(x,y). Now you know the vector potential (which is proportional to the current distribution) from each strip due solely to the current in that strip. Save A2 to disk as well.
Now put both strips in the same space. Solve ∇^2 A = (A + A1 + A2)/λ^2. Here, A is the field you’re solving for, and A1 and A2 are the fixed functions you found before.
I feel like this will work because the London equations are linear, so the solutions to the different problems can just be added together. Of course, when you have the overall solution you can see if it satisfies the London equation.
OK, I realized that the correct boundary condition is not A = 0 on a circular boundary, but rather that the normal derivative on the circular boundary is consistent with that of a current-carrying wire. In other words, use boundary condition
int1d(Th, circle)(mu0I0w/(2piR))
where “circle” is the label of the outer circular boundary and R is the boundary radius.
Thank you very much for your message. The boundary condition that you suggested is given for Bθ (in a system with cylindrical symmetry Br is 0). But as soon as you have another wire (or more) this symmetry is broken.
Initially, I applied the boundary condition
A=+μ0.I/(2π). log(R+)- μ0.I/(2π). log(R-)
where R+ and R- are distances from the running point on the boundary to each center of the couple of wires, one fed with +I and the other fed with -I. The results do not seem realistic.
In addition, I think the linearity of London equations applies to a given system of conductors = if you can solve the problem for (+I, 0) current distribution in a set of 2 conductors, then you can solve (0,-I) with the same method, for the same couple of conductors and apply the superposition theorem to obtain the solution for (+I, -I) current distribution, always for the same set of 2 conductors.
Basically, I am mostly interested in the current distribution inside the superconductors, and marginally in the empty space between the superconductors. And I feel that the current distribution may be sensitive to the boundary conditions that you use. In the physical problem, one applies a total current. I would better trust a solution starting from something equivalent. Intuitively, I expect that the current density is larger on the surface of the superconducting wires (this is the case for a single wire), but with a second wire, the current density is increased for the parts of the surfaces facing the other superconductor and reduced for the parts of the surfaces opposite to the other superconductor. All this is qualitative; I would like to obtain quantitative evaluations.
Thank you again for any help.
Denis
Hi Denis: I think that the way you can solve this problem is as follows. But it would take some tests to be sure.
-
Separately solve for each current-carrying wire, using our agreed-upon BC for a single wire. Save the resulting current distributions J1(x, y), J2(x, y), … for each wire.
-
Now put (say) two wires in the same space. Make the current density in the first wire J1(x, y) and that in the second wire J2(x, y). Then you need to solve
∇^2 A - A/λ^2 = mu_0*(J1 + J2).
Current distributions J1 and J2 are contained in separate FE spaces.
- Now solve the two-wire case with the BC’s equaling the sum of the original BC’s that you used to find J1 and J2. The fact that the two wires are no longer centered in the circular outer space shouldn’t matter much if the outer radius is large.
See also my comments below.
denis.crete
Denis CrétéApril 30Thank you very much for your message. The boundary condition that you suggested is given for Bθ (in a system with cylindrical symmetry Br is 0). But as soon as you have another wire (or more) this symmetry is broken.
Initially, I applied the boundary condition
A=+μ0.I/(2π). log(R+)- μ0.I/(2π). log(R-)
where R+ and R- are distances from the running point on the boundary to each center of the couple of wires, one fed with +I and the other fed with -I. The results do not seem realistic.
I don’t think this BC can work. If R+ and R- are large then the BC will be close to zero, and the amount of current in each wire is only weakly determined by the BC. I think my idea above, to solve for each wire separately, then fix the current in each wire when you put them both in the same space, will work better.
In addition, I think the linearity of London equations applies to a given system of conductors = if you can solve the problem for (+I, 0) current distribution in a set of 2 conductors, then you can solve (0,-I) with the same method, for the same couple of conductors and apply the superposition theorem to obtain the solution for (+I, -I) current distribution, always for the same set of 2 conductors.
Basically, I am mostly interested in the current distribution inside the superconductors, and marginally in the empty space between the superconductors. And I feel that the current distribution may be sensitive to the boundary conditions that you use. In the physical problem, one applies a total current. I would better trust a solution starting from something equivalent. Intuitively, I expect that the current density is larger on the surface of the superconducting wires (this is the case for a single wire), but with a second wire, the current density is increased for the parts of the surfaces facing the other superconductor and reduced for the parts of the surfaces opposite to the other superconductor. All this is qualitative; I would like to obtain quantitative evaluations.
The way to think about this is as follows. Suppose the current distribution in a single current-carrying wire (with no other wire nearby) is J_I(x, y). Now take that wire, with no current in it, and place it in some known magnetic field. The current distribution will be J_B(x, y). Now if you take the wire with the current and also in the feld, the current distribution will be J_I + J_B, by superposition.
If the “external” field actually comes from a fixed current in a second wire, the argument is the same.
Hello,
One of my very first question was, briefly stated “how can I FIX the total current in a conductor” ? My attempts with BC gave results highly dependent on the domain shape and essentially not reliable. I did not succeed either in adding a Lagrange factor to the stiffness matrix. Thank you if you have a method.
A totally agree with you for applying the superposition theorem with/without external B and with/without applied current to a single wire. But if the “external” magnetic field comes from a nearby second (superconducting) wire, the current distribution in both superconducting wires is changed (and this is precisely this change that I would like to evaluate…); in the approach you suggest, I understand that you impose the current distribution in the last step, i.e. when we add the second wire.
I have a few other questions :
As in a superconductor A= - λ2μ0×J, should’nt the equation be written either:
∇2 A= - μ0×(J1 + J2)
or
∇2 A - A/λ2 =0
?
Thus I don’t understand why solving ∇2 A - A/λ2 = - μ0×(J1 + J2) will give you the solution of ∇2 A - A/λ2 =0 .
As the currents I am running in each wire are opposite, the BC that I would use for the single-wire case are opposite. Then, at step 3, the sum of the BC is identically zero if you neglect the fact that the wires are not perfectly centered.
Writing the BC as : A=+μ0.I/(2π). log(R+)- μ0.I/(2π). log(R-)
was an attempt to make the sum of the BC, taking into account the fact that the wires were not at the same location. And the result is not reliable.
Please tell me if I misunderstood your suggestions. Thank you anyway for trying to explain.
Best regards
Denis
Hi Denis: I am starting to think that my method of superposition won’t work. I did, however, find a paper in the literature that addresses exactly the problem you’re trying to solve. It’s
https://www-sciencedirect-com.ezproxy2.library.colostate.edu/science/article/pii/0021999179900706
I’m reading it carefully now and trying to understand some of the subtleties, especially concerning their function X. I think it may be related to the fact that for infinitely long wires, any flux between the wires needs to satisfy the fluxoid condition. I tried a simple simulation with two wires and the BCs they discussed (but without the function X), and it didn’t work. So now I will try with the X-function, although this may require an iterative solution, since X depends on the integral of A over the surface of each wire.
Hello Stuart,
Thank you very much for the link; although I do not have access to the library, I could reformulate the search on the net and find 2 interesting papers which make analytical evaluations for the current density in a system of superconducting wires:
D. Oswald, “Current and current density distribution in parallel superconducting wires in a planar array in the Meissner state”, CRYOGENICS. MAY 1973, p290
and
A.R. Sass et al., “CURRENT DISTRIBUTION IN A THIN FILM SUPERCONDUCTING STRIP TRANSMISSION LINE”, NASA Technical note D-3344, march 1966 https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19660010306.pdf
I think we can calculate this current distribution with the method used in these publications. But probably with little help from FreeFEM.
I could obtain some results with InductEx v4.23 (or 4.30, I am not sure), but there is a limitation on the complexity of the circuit. I expect that the licenced version of InductEx (5.xx) can model the 3D circuit. This software uses Fasthenry to evaluate the inductances of (superconducting) circuits.
Could you please send me the reference of the paper you were mentioning? I am still interested in learning about the approach.
And still a bit frustrated to feel that I can’t use FreeFEM correctly…
Best regards
Denis
In a private message you asked me the question below. I include it here in case anyone in the community is interested. Your question is about the utility or meaning of the variable X. I will try to answer this as best I can. The gauge-independent version of the London equation is (in Alsop’s units; see Eq. 2.9) ∇ x J = –B. Because B = ∇ x A, we can write London’s equation as ∇ x (J + A) = 0. Now if a field (such as J + A) has zero curl, then it mist be true that the field can be written as the gradient of a scalar potential, J + A = -∇ϕ (this is Eq. 2.11). If you take the divergence of both sides of this equation you get ∇·J + ∇·A = -∇^2ϕ. But ∇·J = 0 by charge conservation, so this tells us that ∇·A = -∇^2ϕ.
In 2D, ∇·A = 0 so that ∇^2ϕ = 0. Thus, as Alsop says, ∇ϕ = X = constant. This leads to Eq. 2.12, ∇^2A - A = X. But notice that this only holds inside the superconductor. Outside we have ∇^2A = 0. So X is a constant inside the superconductor, but zero outside. So in terms of the FEM solution, it is not an overall additive constant, but rather a function of position. Thus it will affect the solution.
I have attempted to solve a problem like the ones in Alsop using FreeFEM. I attach it to this reply. I set up two circular wires and put current I0 through one. However, I can’t get the current through the other (which does not have a current enforced in it) to be zero. I would think that if you had two parallel wires, one with a current but the other floating, the second wire would tend not to have a current. But this may be a peculiarity of a 2D simulation. Infinitely long wires with current have an infinite energy, for example. I note that in Alsop, all they geometries have a ground plane, which gets rid of this issue. However, I tried a FreeFEM simulation of two wires with aground plane, and the second wire still had a current in it.
You can try out this edp file and see what you think.
Your earlier message:
"Thank you very much for the paper and for the time you spent to help me !
I looks indeed very close to what I am trying to do with FreeFEM.
What I do not understand in this paper, is the utility of variable X. This variable is the degree of freedom which is still left after choosing the Coulomb gauge (∇.A=0);
and it is said that it is a constant in space.
Besides the fact that London equation, established for the Coulomb gauge, states A=-μ0 λ² j (i.e. in eq. 2.11, X=0…), solving ΔA - A/λ² = X (eq. 2.12), is equivalent to solving ΔA - A/λ²= 0, where A=A+λ²X… It just changes the value of the vector potential by a constant… This should not change anything in the calculations… provided they comply with “constant” behaviour of X.
It becomes very strange in eq.14 where X1…Xn are introduced for the set of n superconducting wires, as if they could be chosen independently. In my understanding, they are linked by the relations (λi)²Xi=constant (i.e. independant of i).
Best regards
Denis"
TwoWires.edp (3.2 KB)
Hello,
Thank you very much! I think you helped me a lot, because after introducing an X also for the second wire (which has no current), I think we obtain the right result. All the Xk must be evaluated at each iteration, as all the superconducting wires are immersed in a space where the potential vector A is non zero (as soon as current flows in at least one wire).
Attached is the code containing the modifications necessary after renaming the “X” parameter for the left wire “XL” and introducing XR for the right wire; plus a few other “cosmetic” changes to check the credibility of the results (I did not try yet to reproduce the results of Alsop).
Just to have more insight on the X parameter, I went back to the equation of quantum transport for spin-0 particles (Cooper pairs are such 0 spin particles) which reads (in SI units)
A+μ0λ2J = Φ0/2π . ∇φ
where
Φ0 is the flux quantum and φ is the phase of the macroscopic wave function of the superconducting state. This equation is gauge-independant.
Therefore,
λ2X= Φ0/2π. ∇φ
In this case X is constant and ∇φ = 2π/Λ (Λ is the wavelength of the superconducting wave function), leading:
λ2X= Φ0/Λ
Note that X and Λ are gauge-dependant. Because of Λ, this equation clearly shows that X is defined separately for each superconducting wire.
Thank you again !
TwoWires-WRMCJB.edp (4.2 KB)