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 = function.stack(['r cos(θ)', 'r sin(θ)', 'z'] @ ns)
 86    ns.x = ns.rz @ ns.rot
 87    ns.define_for('x', gradient='∇', jacobians=('dV', 'dS'), curl='curl')
 88
 89    ns.basis = RZ.basis('spline', degree=degree, removedofs=[[0, -1], [-1]])
 90    ns.A = function.dotarg('A', ns.basis, dtype=complex) * ns.eθ
 91    ns.Atest = function.dotarg('Atest', ns.basis, dtype=float) * ns.eθ
 92    ns.B_i = 'curl_ij(A_j)'
 93    ns.E_i = '-j ω A_i'
 94    ns.Jind_i = 'σ E_i'
 95    ns.I = 1
 96    ns.Jext_i = 'eθ_i I / π rwire^2'
 97    ns.J_i = 'Jext_i + Jind_i'
 98
 99    res = REV.integral(RZ.integral('-∇_j(Atest_i) ∇_j(A_i) dV' @ ns, degree=2*degree), degree=0)
100    res += REV.integral(RZ['coil'].integral('μ0 Atest_i J_i dV' @ ns, degree=2*degree), degree=0)
101
102    state = solver.solve_linear(('A',), (res.derivative('Atest'),))

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.

108    ns.Borig = ns.B
109    ns.Bbasis = RZ.basis('spline', degree=degree)
110    ns.B = function.dotarg('B', ns.Bbasis, dtype=complex, shape=(2,)) @ ns.rot
111    ns.Btest = function.dotarg('Btest', ns.Bbasis, shape=(2,)) @ ns.rot
112    res = REV.integral(RZ.integral('Btest_i (B_i - Borig_i) dV' @ ns, degree=2*degree), degree=0)
113    state = solver.solve_linear(('B',), (res.derivative('Btest'),), arguments=state)
114
115    with export.mplfigure('magnetic-potential-1.png', dpi=300) as fig:
116        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.

119        smpl = REV0 * RZ[:-1, :-1].sample('bezier', 5)
120        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, **state)
121        Amax = abs(A).max()
122        Bmax = abs(Bmag).max()
123        levels = numpy.linspace(-Amax, Amax, 32)[1:-1]
124        r = numpy.concatenate([r, r], axis=0)
125        z = numpy.concatenate([z, -z], axis=0)
126        A = numpy.concatenate([A, A], axis=0)
127        Bmag = numpy.concatenate([Bmag, Bmag], axis=0)
128        tri = numpy.concatenate([smpl.tri+i*smpl.npoints for i in range(2)])
129        hull = numpy.concatenate([smpl.hull+i*smpl.npoints for i in range(2)])
130        imBi = ax.tripcolor(-r, z, tri, Bmag.imag, shading='gouraud', cmap='Greens')
131        imBi.set_clim(0, Bmax)
132        ax.tricontour(-r, z, tri, -A.imag, colors='k', linewidths=.5, levels=levels)
133        imBr = ax.tripcolor(r, z, tri, Bmag.real, shading='gouraud', cmap='Greens')
134        imBr.set_clim(0, Bmax)
135        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.

138        smpl = REV0 * RZ['coil'].sample('bezier', 5)
139        r, z, J = smpl.eval(['r', 'z', 'J_1'] @ ns, **state)
140        Jmax = abs(J).max()
141        r = numpy.concatenate([r, r], axis=0)
142        z = numpy.concatenate([z, -z], axis=0)
143        J = numpy.concatenate([J, J], axis=0)
144        tri = numpy.concatenate([smpl.tri+i*smpl.npoints for i in range(2)])
145        imJi = ax.tripcolor(-r, z, tri, -J.imag, shading='gouraud', cmap='bwr')
146        imJi.set_clim(-Jmax, Jmax)
147        imJr = ax.tripcolor(r, z, tri, J.real, shading='gouraud', cmap='bwr')
148        imJr.set_clim(-Jmax, Jmax)

Minor ticks at element edges.

150        rticks = RZ.boundary['bottom'].interfaces.sample('gauss', 0).eval(ns.r)
151        ax.set_xticks(numpy.concatenate([-rticks[::-1], [0], rticks]), minor=True)
152        zticks = RZ.boundary['left'].interfaces.sample('gauss', 0).eval(ns.z)
153        ax.set_yticks(numpy.concatenate([-zticks[::-1], [0], zticks]), minor=True)
154        ax.tick_params(direction='in', which='minor', bottom=True, top=True, left=True, right=True)

Real and imag indicator.

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

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

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

Minor ticks at element edges.

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