coil.py

In this script we compute the magnetic field induced by a DC or AC current in one or several toroidal conductors. This problem is modeled with the quasi-static magnetic vector potential with Lorenz gauge:

\[∇_j(∇_j(A_i)) = -μ_0 J_i,\]

where \(A\) is the magnetic vector potential, \(J\) the current density and \(μ_0\) the magnetic permeability. The magnetic field \(B\) is then given by the curl of the magnetic vector potential. The current density is the sum of an external current \(J_\text{ext}\) and the current induced by the magnetic field, \(J_\text{ind}\). The external current is given by

\[J_{\text{ext};i} = \frac{I}{π r_\text{wire}^2} \cos(ω t) e_{θi},\]

inside the conductor and zero everywhere else, where \(ω = 2 π f\). The induced current follows from Faraday’s law of induction and Ohm’s law:

\[J_{\text{ind};i} = -σ ∂_t(A_i),\]

where \(σ\) is the conductivity, which is non-zero only inside the conductor.

We can solve the temporal component of \(A\) by letting \(A_i = \Re(\hat{A}_i \exp{j ω t})\). This problem in terms of \(\hat{A}\) is:

\[∇_j(∇_j(\hat{A}_i) = -μ_0 \hat{J}_i,\]

with

\[\hat{J}_{\text{ext};i} = \frac{I}{π r_\text{wire}^2} e_{θi},\]

and

\[\hat{J}_{\text{ind};i} = -j ω σ \hat{A}_i.\]
43from nutils import cli, export, function, mesh, solver, testing
44from nutils.expression_v2 import Namespace
45import functools
46import numpy
47
48
49def main(nelems: int = 50, degree: int = 3, freq: float = 0., nturns: int = 1, rwire: float = .0025, rcoil: float = 0.025):
50
51    ns = Namespace()
52    ns.j = 1j
53    ns.π = numpy.pi
54    ns.f = freq
55    ns.ω = '2 π f'
56    ns.μ0 = '4e-7 π'  # magnetic permeability in vacuum
57    ns.σ = 5.988e7  # conductivity of copper
58    ns.rcoil = rcoil  # radius of coil (to the center of the wire)
59    ns.rwire = rwire  # radius of the wire

The problem is axisymmetric in the z-axis and symmetric in z=0. We create a 2D domain RZ covering the quarter plane [0,inf)^2, multiply this with a 1D domain REV with one element spanning [-π,π] and transform the geometry and vector bases from cylindrical to cartesian coordinates. The RZ domain is mapped from [0,1] to [0,inf) using an arctanh transform. Finally, a Neumann boundary condition is used at z=0 to obtain symmetry in z=0.

68    RZ, ns.rz0 = mesh.rectilinear([numpy.linspace(0, 1, nelems)]*2, space='RZ')
69    REV, ns.θ = mesh.line([-numpy.pi, numpy.pi], bnames=['start', 'end'], space='Θ')
70    REV0 = REV.refined[:1].boundary['end'].sample('bezier', 2)
71    ns.rz = numpy.arctanh(ns.rz0) * 2 * rcoil
72    ns.r, ns.z = ns.rz

Trimming of the wires. The centers of the wires are located at rcoil,zwires.

77    ns.zwires = (numpy.arange(nturns) - (nturns - 1) / 2) * 4 * rwire
78    ns.dwires = ns.rwire - numpy.sqrt((ns.r - ns.rcoil)**2 + functools.reduce(function.min, (ns.z - ns.zwires)**2))
79    RZ = RZ.withsubdomain(coil=RZ[:-1, :-1].trim(ns.dwires/ns.rwire, maxrefine=4))
80
81    ns.rot = function.stack([function.scatter(function.trignormal(ns.θ), 3, [0, 1]), function.kronecker(1., 0, 3, 2)])
82    ns.eθ = function.stack(['-sin(θ)', 'cos(θ)', '0'] @ ns)
83
84    X = RZ * REV
85    ns.x = ns.rz @ ns.rot
86    ns.define_for('x', gradient='∇', jacobians=('dV', 'dS'), curl='curl')
87    ns.add_field(('A', 'Atest'), RZ.basis('spline', degree=degree, removedofs=[[0, -1], [-1]])[:,numpy.newaxis] * ns.eθ, dtype=complex)
88    ns.B_i = 'curl_ij(A_j)'
89    ns.E_i = '-j ω A_i'
90    ns.Jind_i = 'σ E_i'
91    ns.I = 1
92    ns.Jext_i = 'eθ_i I / π rwire^2'
93    ns.J_i = 'Jext_i + Jind_i'
94
95    res = REV.integral(RZ.integral('-∇_j(Atest_i) ∇_j(A_i) dV' @ ns, degree=2*degree), degree=0)
96    res += REV.integral(RZ['coil'].integral('μ0 Atest_i J_i dV' @ ns, degree=2*degree), degree=0)
97
98    args = solver.solve_linear('A:Atest,', res)

Since the coordinate transformation is singular at r=0 we can’t evaluate B (the curl of A) at r=0. We circumvent this problem by projecting B on a basis.

104    ns.Borig = ns.B
105    ns.add_field(('B', 'Btest'), RZ.basis('spline', degree=degree), ns.rot, dtype=complex)
106    res = REV.integral(RZ.integral('Btest_i (B_i - Borig_i) dV' @ ns, degree=2*degree), degree=0)
107    args = solver.solve_linear('B:Btest,', res, arguments=args)
108
109    with export.mplfigure('magnetic-potential-1.png', dpi=300) as fig:
110        ax = fig.add_subplot(111, aspect='equal', xlabel='$x_0$', ylabel='$x_2$', adjustable='datalim')

Magnetic vector potential. r < 0 is the imaginary part, r > 0 the real part.

113        smpl = REV0 * RZ[:-1, :-1].sample('bezier', 5)
114        r, z, A, Bmag = smpl.eval(['r', 'z', 'A_1', 'sqrt(real(B_i) real(B_i)) + sqrt(imag(B_i) imag(B_i)) j'] @ ns, **args)
115        Amax = abs(A).max()
116        Bmax = abs(Bmag).max()
117        levels = numpy.linspace(-Amax, Amax, 32)[1:-1]
118        r = numpy.concatenate([r, r], axis=0)
119        z = numpy.concatenate([z, -z], axis=0)
120        A = numpy.concatenate([A, A], axis=0)
121        Bmag = numpy.concatenate([Bmag, Bmag], axis=0)
122        tri = numpy.concatenate([smpl.tri+i*smpl.npoints for i in range(2)])
123        hull = numpy.concatenate([smpl.hull+i*smpl.npoints for i in range(2)])
124        imBi = ax.tripcolor(-r, z, tri, Bmag.imag, shading='gouraud', cmap='Greens')
125        imBi.set_clim(0, Bmax)
126        ax.tricontour(-r, z, tri, -A.imag, colors='k', linewidths=.5, levels=levels)
127        imBr = ax.tripcolor(r, z, tri, Bmag.real, shading='gouraud', cmap='Greens')
128        imBr.set_clim(0, Bmax)
129        ax.tricontour(r, z, tri, A.real, colors='k', linewidths=.5, levels=levels)

Current density (wires only). r < 0 is the imaginary part, r > 0 the real part.

132        smpl = REV0 * RZ['coil'].sample('bezier', 5)
133        r, z, J = smpl.eval(['r', 'z', 'J_1'] @ ns, **args)
134        Jmax = abs(J).max()
135        r = numpy.concatenate([r, r], axis=0)
136        z = numpy.concatenate([z, -z], axis=0)
137        J = numpy.concatenate([J, J], axis=0)
138        tri = numpy.concatenate([smpl.tri+i*smpl.npoints for i in range(2)])
139        imJi = ax.tripcolor(-r, z, tri, -J.imag, shading='gouraud', cmap='bwr')
140        imJi.set_clim(-Jmax, Jmax)
141        imJr = ax.tripcolor(r, z, tri, J.real, shading='gouraud', cmap='bwr')
142        imJr.set_clim(-Jmax, Jmax)

Minor ticks at element edges.

144        rticks = RZ.boundary['bottom'].interfaces.sample('gauss', 0).eval(ns.r)
145        ax.set_xticks(numpy.concatenate([-rticks[::-1], [0], rticks]), minor=True)
146        zticks = RZ.boundary['left'].interfaces.sample('gauss', 0).eval(ns.z)
147        ax.set_yticks(numpy.concatenate([-zticks[::-1], [0], zticks]), minor=True)
148        ax.tick_params(direction='in', which='minor', bottom=True, top=True, left=True, right=True)

Real and imag indicator.

150        spine = next(iter(ax.spines.values()))
151        ax.axvline(0, color='k')
152        ax.text(0, .95, '← imag ', horizontalalignment='right', verticalalignment='bottom', transform=ax.get_xaxis_transform())
153        ax.text(0, .95, ' real →', horizontalalignment='left', verticalalignment='bottom', transform=ax.get_xaxis_transform())
154        ax.set_xlim(-2*rcoil, 2*rcoil)
155        ax.set_ylim(-2*rcoil, 2*rcoil)
156        fig.colorbar(imJr, label='$J_1$')
157        fig.colorbar(imBr, label='$|B|$')
158
159    if freq == 0:
160        ns.δ = function.eye(3)

Reference: https://physics.stackexchange.com/a/355183

162        ns.Bexact = ns.δ[2] * ns.μ0 * ns.I * ns.rcoil**2 / 2 * ((ns.rcoil**2 + (ns.z - ns.zwires)**2)**(-3/2)).sum()
163        smpl = REV0 * RZ[:-1, :-1].boundary['left'].sample('bezier', 5)
164        B, Bexact, z = smpl.eval(['real(B_2)', 'Bexact_2', 'z'] @ ns, **args)
165        z = numpy.concatenate([-z[::-1], z])
166        B = numpy.concatenate([B[::-1], B])
167        Bexact = numpy.concatenate([Bexact[::-1], Bexact])
168        with export.mplfigure('magnetic-field-x2-axis.png', dpi=300) as fig:
169            ax = fig.add_subplot(111, xlabel='$x_2$', ylabel='$B_2$', title='$B_2$ at $x_0 = x_1 = 0$')
170            ax.plot(z, B, label='FEM')
171            ax.plot(z, Bexact, label='exact', linestyle='dotted')
172            ax.legend()

Minor ticks at element edges.

174            zticks = RZ.boundary['left'].interfaces.sample('gauss', 0).eval(ns.z)
175            ax.set_xticks(numpy.concatenate([-zticks[::-1], [0], zticks]), minor=True)
176            ax.tick_params(axis='x', direction='in', which='minor', bottom=True, top=True)
177            ax.set_xlim(-2*rcoil, 2*rcoil)
178
179    return args
180
181
182class test(testing.TestCase):
183
184    def test_dc(self):
185        args = main(nelems=16, degree=2)
186        with self.subTest('A.real'):
187            self.assertAlmostEqual64(args['A'].real, '''
188                eNoNke9rzWEYh5NzVmtnvud5nvv+3PdzTn7lIIRlL3Rq/wArinFGaytFo6xjTedISMwsJsNksbJYtlIS
189                U9pqLcqJKL9ytL3xYm92kpkQ2vL9B67P9el6TS/oHuVpPb13zW7WZu2U2WaG4t8CF8xWVsS+YgZF3MYu
190                /OYLTHyFyijrXllrNxvEqxaVa1S/yJBk5CfaEUMnz1MzPXcxV23JVAWjOq4D2qAL9YakZBAp9HKE99F9
191                99E+NcWgw5/yaT+tJzWm3WLlEiI4wu9oKdW6TTYTL/m//oPf4T9rvU7IXvmE7RjjFB+lAXfZjsRrk2uT
192                qxM3fcSfDTfaJSqn8YubeJhKbtIG5kdiImESHX5ez2iFXpWk9MHjPE/Rckq4jDnhO/0xv8SPhfwZOScq
193                d7EG/VzGW0ODHvNdS+GDa7pTy/WJNMgcrmMlBln4ALW6l2aZrtCk/pO3cksaRaSALCrRx8pt1OX+mLzk
194                5JDUS01ILmEYOWzEJB/nKGep1y22j/AYD3AH3chjD6oRxRu+yDVcpN3U49K2wAV+xqP8kPu5i9u4jjfw
195                HI1Ta9ihya2zLdRCh+kg7adGqqMtlKZVFKNpN+JyboFL2f8Z6oV2''')
196
197    def test_ac_5(self):
198        args = main(nelems=16, degree=2, freq=1000, nturns=5)
199        with self.subTest('A.imag'):
200            self.assertAlmostEqual64(args['A'].imag, '''
201                eNoNkEtIlGEYhRcWBqVWNgsxujBImxJmIIwEc2FJCYVaQxEErRIvkyFOhqTB4CLSpE3TxQiLCjXQEdJN
202                QUkZKMxCyxYhmpvvPe/7ft//K11AvPSvzuJZnPOcOZ3TeV1SX3PtERu3n2zEfXRNXqkfXU6s9P5O8wiP
203                8nue5kXeLhXyRHL0mVbZApfn1fj1K4MYwgjeYQLzWEcZpzhPXkqtHrAhF/Pqlzdpg9YpG/kIowb3wGjg
204                rTImSU3YUTfhN9MNaqcEdVIfDdIaNeAPUnxGQpplS12Fn0+7KItCVBhkLSVpC17jNG/wFxnWRbvgtZmk
205                +WxumQ4zY6rNX/ODmrGfp7hH4vrIdnuPzQuTMU9Nv/lpeswe+h407Aic6qRcr9iT3jE6SoepjE5RJXXR
206                MOXiPkI8xBfkoEbtNu859dNAsGycvhFRGHFM4Xjwx3nZqbvtrCtCEQ4hivLALoE+zGIvt/IvvhbwTX3j
207                0khjDB8wjQX8QyFXcidP8j7plbCuaZeLcYwv8VVu5HZ+wG85w6sckZsyI+c0xza6YimWiJTICamSy3Jd
208                7spAwLL1rKa12t5xLdqirdqmtzWp3ZrSVzquGXVaYC/ar/ah+w/zsU82''')
209
210
211if __name__ == '__main__':
212    cli.run(main)