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:
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
inside the conductor and zero everywhere else, where \(ω = 2 π f\). The induced current follows from Faraday’s law of induction and Ohm’s law:
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:
with
and
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)