From 7842d0545e776702a6cea05de4bf6a14cd2707bf Mon Sep 17 00:00:00 2001 From: simbilod Date: Wed, 12 Jun 2024 23:39:15 -0700 Subject: [PATCH 1/4] add refs --- docs/references.bib | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/docs/references.bib b/docs/references.bib index b2c88001..8088b875 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -419,3 +419,29 @@ @article{Williams1997 month = jul, pages = {405} } +@book{Haus1989, + title = {Electromagnetic Fields and Energy}, + author = {Haus, H.A. and Melcher, J.R.}, + isbn = {9780132492775}, + lccn = {88007281}, + url = {https://ocw.mit.edu/courses/res-6-001-electromagnetic-fields-and-energy-spring-2008/}, + year = {1989}, + publisher = {Prentice Hall} +} +@book{Larson2013, + title = {The {{Finite Element Method}}: {{Theory}}, {{Implementation}}, and {{Applications}}}, + shorttitle = {The {{Finite Element Method}}}, + author = {Larson, Mats G. and Bengzon, Fredrik}, + year = {2013}, + series = {Texts in {{Computational Science}} and {{Engineering}}}, + volume = {10}, + publisher = {Springer}, + address = {Berlin, Heidelberg}, + doi = {10.1007/978-3-642-33287-6}, + urldate = {2024-06-12}, + copyright = {https://www.springernature.com/gp/researchers/text-and-data-mining}, + isbn = {978-3-642-33286-9 978-3-642-33287-6}, + langid = {english}, + keywords = {adaptive finite element methods,applications of finite elements,error estimates,finite element method,implementation of finite elements,partial differential equations}, + file = {C:\Users\simon\Zotero\storage\AQGUB2BJ\Larson and Bengzon - 2013 - The Finite Element Method Theory, Implementation,.pdf} +} From c5b59386bb60811367ca1e7ce0080b819eb589c1 Mon Sep 17 00:00:00 2001 From: simbilod Date: Wed, 12 Jun 2024 23:39:55 -0700 Subject: [PATCH 2/4] inductor example --- docs/electronics/examples/inductor.py | 378 ++++++++++++++++++++++++++ femwell/magnetostatic.py | 207 ++++++++++++++ 2 files changed, 585 insertions(+) create mode 100644 docs/electronics/examples/inductor.py create mode 100644 femwell/magnetostatic.py diff --git a/docs/electronics/examples/inductor.py b/docs/electronics/examples/inductor.py new file mode 100644 index 00000000..cb305fe4 --- /dev/null +++ b/docs/electronics/examples/inductor.py @@ -0,0 +1,378 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.16.2 +# kernelspec: +# display_name: waveguidesmodes +# language: python +# name: python3 +# --- + +# # Magnetostatics and inductance +# +# We follow the theory outlined in the documentation to compute the magnetic field of various currents distributions. + +# + +import matplotlib.pyplot as plt +import numpy as np +import shapely +from meshwell.model import Model +from meshwell.polysurface import PolySurface +from scipy.interpolate import griddata +from skfem import ( + Basis, + BilinearForm, + ElementDG, + ElementTriN1, + ElementTriP0, + ElementTriP1, + ElementVector, + FacetBasis, + Functional, + LinearForm, + asm, + condense, + solve, +) +from skfem.helpers import curl, dot, grad, inner +from skfem.io import from_meshio +from skfem.io.meshio import from_meshio +from skfem.models.general import curluv +from tqdm import tqdm + +from femwell.magnetostatic import solve_magnetostatic_2D +from femwell.mesh import mesh_from_OrderedDict +from femwell.visualization import plot_domains, plot_subdomain_boundaries + +# %load_ext autoreload +# %autoreload 2 + + +# + +core_permeability = 1 +diameter = 2 +length = 5 +current = 1 + +dr_domain = 50 +dz_domain = 50 +# - + + +# Make a mesh +# +#
+# Note: below we use meshwell instead of femwell's built-in backend. You can install meshwell with `pip install meshwell` +#
+ + +# + tags=["remove-stderr", "hide-output"] +def continuous_solenoid_cross_section_mesh( + diameter=diameter, + length=length, + dr_domain=10, + dz_domain=10, +): + + solenoid_polygon = shapely.box(-diameter / 2, -length / 2, diameter / 2, length / 2) + left_air_polygon = shapely.box(-dr_domain / 2, -dz_domain / 2, -diameter / 2, dz_domain / 2) + right_air_polygon = shapely.box(diameter / 2, -dz_domain / 2, dr_domain / 2, dz_domain / 2) + center_air_polygon = shapely.box(-diameter / 2, -dz_domain / 2, diameter / 2, dz_domain / 2) + + model = Model() + + solenoid = PolySurface( + polygons=solenoid_polygon, + model=model, + physical_name="solenoid", + mesh_bool=True, + resolution={"resolution": 0.1, "DistMax": 2, "SizeMax": 0.2}, + mesh_order=1, + ) + left_air = PolySurface( + polygons=left_air_polygon, + model=model, + physical_name="left", + mesh_bool=True, + mesh_order=2, + ) + right_air = PolySurface( + polygons=right_air_polygon, + model=model, + physical_name="right", + mesh_bool=True, + mesh_order=2, + ) + center_air = PolySurface( + polygons=center_air_polygon, + model=model, + physical_name="center", + mesh_bool=True, + mesh_order=2, + ) + + return from_meshio( + model.mesh( + entities_list=[solenoid, left_air, center_air, right_air], + filename="mesh.msh", + default_characteristic_length=0.5, + ) + ) + + +mesh = continuous_solenoid_cross_section_mesh( + diameter=diameter, + length=length, + dr_domain=dr_domain, + dz_domain=dz_domain, +) +mesh.draw().show() + +plot_domains(mesh) +plt.show() + +plot_subdomain_boundaries(mesh) +# - + + +# We use these three background volumes to easily tag the left and right sides of the solenoid to define currents + + +# + +def magnetic_potential(mesh, current=1, permeability_mu_r=1): + + # Permeability defined constant-wise to handle sharp discontinuities + basis_mu_r = Basis(mesh, ElementTriP0()) + mu_r = basis_mu_r.ones() + mu_r[basis_mu_r.get_dofs(elements=("solenoid"))] = permeability_mu_r + + basis_A, A = solve_magnetostatic_2D( + basis_mu_r=basis_mu_r, + mu_r=mu_r, + mesh=mesh, + current_densities={ + "solenoid___left": -current / length / 2, # convert current to current density + "solenoid___right": current / length / 2, + }, + ) + + return basis_A, A, basis_mu_r, mu_r + + +basis_A, A, basis_mu_r, mu_r = magnetic_potential( + mesh, current=current, permeability_mu_r=core_permeability +) +# - + +# The magnetic field is $\vec{B} = \nabla \times \vec{A}$, or in terms of test functions, + +# + +vector_basis = Basis(mesh, ElementVector(ElementTriP1())) + +# Need a minus sign here? +bfield = -1 * asm( + LinearForm(curluv).partial(basis_A.interpolate(A)), + vector_basis, +) + +# + +fig, ax = plt.subplots() +for subdomain in ["solenoid"]: + basis_mu_r.mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True) +basis_A.plot(A, ax=ax, shading="gouraud", colorbar=True) + +# Define the grid for streamplot interpolation +Nx, Ny = 100, 100 +xmin, xmax = mesh.p[0].min(), mesh.p[0].max() +ymin, ymax = mesh.p[1].min(), mesh.p[1].max() +grid_x, grid_y = np.meshgrid(np.linspace(xmin, xmax, Nx), np.linspace(ymin, ymax, Ny)) + +# Interpolate the magnetic field on the grid +bfield_x = bfield.reshape((-1, 2))[:, 0] +bfield_y = bfield.reshape((-1, 2))[:, 1] +grid_bfield_x = griddata(mesh.p.T, bfield_x, (grid_x, grid_y), method="cubic") +grid_bfield_y = griddata(mesh.p.T, bfield_y, (grid_x, grid_y), method="cubic") + +# Plot the streamlines +lw = ( + 5 + * np.sqrt(grid_bfield_x**2 + grid_bfield_y**2) + / np.max(np.sqrt(grid_bfield_x**2 + grid_bfield_y**2)) +) +ax.streamplot(grid_x, grid_y, grid_bfield_x, grid_bfield_y, linewidth=lw, density=1.2, color="k") +plt.show() +# - + +max_A_femwell = np.max(A) + +min_A_femwell = np.min(A) + +max_B_femwell = np.max(np.sqrt(grid_bfield_x**2 + grid_bfield_y**2)) + +min_B_femwell = np.min(np.sqrt(grid_bfield_x**2 + grid_bfield_y**2)) + +# We can compare this to [the analytical solution](https://en.wikipedia.org/wiki/Solenoid#Finite_continuous_solenoid): + +# + +from scipy.special import ellipe, ellipk, elliprf, elliprj + + +# See https://github.com/scipy/scipy/issues/4452 +def ellipp(n, m): + if m >= 1: + raise ValueError("m must be < 1") + y = 1 - m + rf = elliprf(0, y, 1) + rj = elliprj(0, y, 1, 1 - n) + return rf + rj * n / 3 + + +def ellippinc(phi, n, m): + nc = np.floor(phi / np.pi + 0.5) + phi -= nc * np.pi + sin_phi = np.sin(phi) + sin2_phi = sin_phi * sin_phi + sin3_phi = sin2_phi * sin_phi + x = 1 - sin2_phi + y = 1 - m * sin2_phi + rf = elliprf(x, y, 1) + rj = elliprj(x, y, 1, 1 - n * sin2_phi) + val = sin_phi * rf + sin3_phi * rj * n / 3 + if nc != 0: + rp = ellipp(n, m) + val += 2 * nc * rp + return val + + +def m(rho, R, eta): + return 4 * R * rho / ((R + rho) ** 2 + eta**2) + + +def n(rho, R): + return 4 * R * rho / (R + rho) ** 2 + + +def continuous_solenoid_analytical_A(rho, z, R, l): + """For mu0 = 1, I = 1""" + + def term(eta): + sqrt_term = np.sqrt((R + rho) ** 2 + eta**2) + m_val = m(rho, R, eta) + n_val = n(rho, R) + k_term = ellipk(m_val) + e_term = ellipe(m_val) + inc_term = np.vectorize(ellippinc)(np.pi / 2, n_val, m_val) + + return ( + eta + / sqrt_term + * ( + (m_val + n_val - m_val * n_val) / (m_val * n_val) * k_term + - 1 / m_val * e_term + + (n_val - 1) / n_val * inc_term + ) + ) + + eta_plus = z + l / 2 + eta_minus = z - l / 2 + + return R / (np.pi * l) * (term(eta_plus) - term(eta_minus)) + + +""" +Conversions + +R = a +h2 = n +k2 = m +""" + + +def continuous_solenoid_analytical_Bz(rho, z, R, l): + """For mu0 = 1, I = 1""" + + def term(eta): + m_val = m(rho, R, eta) + n_val = n(rho, R) + k_term = ellipk(m_val) + inc_term = np.vectorize(ellippinc)(np.pi / 2, n_val, m_val) + + return eta * np.sqrt(m_val) * (k_term + (R - rho) / (R + rho) * inc_term) + + eta_plus = z + l / 2 + eta_minus = z - l / 2 + + return 1 / (4 * np.pi * l * np.sqrt(R * rho)) * (term(eta_plus) - term(eta_minus)) + + +def continuous_solenoid_analytical_Br(rho, z, R, l): + """For mu0 = 1, I = 1""" + + def term(eta): + m_val = m(rho, R, eta) + k_term = ellipk(m_val) + e_term = ellipe(m_val) + + return (m_val - 2) / np.sqrt(m_val) * k_term + 2 / np.sqrt(m_val) * e_term + + eta_plus = z + l / 2 + eta_minus = z - l / 2 + + return np.sqrt(R) / (2 * np.pi * l * np.sqrt(rho)) * (term(eta_plus) - term(eta_minus)) + + +zs = np.linspace(-dz_domain / 2, dz_domain / 2, 200) +rhos = np.linspace(-dr_domain / 2, dr_domain / 2, 200) +rhos_positive = rhos[rhos >= 0] + +X, Y = np.meshgrid(rhos, zs) +X_positive, Y_positive = np.meshgrid(rhos_positive, zs) +A = continuous_solenoid_analytical_A(X, Y, R=diameter / 2, l=length) + +Bz_positive = continuous_solenoid_analytical_Bz(X_positive, Y_positive, R=diameter / 2, l=length) +Br_positive = continuous_solenoid_analytical_Br(X_positive, Y_positive, R=diameter / 2, l=length) + +fig, ax = plt.subplots() +stream = ax.streamplot( + X_positive, Y_positive, Br_positive, Bz_positive, color="black", linewidth=0.5, density=0.75 +) +ax.set_xlabel("rho") +ax.set_ylabel("z") +ax.set_title("Streamplot of Magnetic Field (Br, Bz) in Continuous Solenoid") + +c = ax.contourf(X, Y, A, levels=50, cmap="jet") +fig.colorbar(c, ax=ax) +ax.set_xlabel("r") +ax.set_ylabel("z") +ax.set_title(r"Analytical $A_{\phi}$ of Continuous Solenoid") +plt.show() +# - + +min_A_analytical = np.min(A) +max_A_analytical = np.max(A) +min_B_analytical = np.min(np.sqrt(Br_positive**2 + Bz_positive**2)) +max_B_analytical = np.max(np.sqrt(Br_positive**2 + Bz_positive**2)) + +min_A_analytical, min_A_femwell + +max_A_analytical, max_A_femwell + +max_B_analytical, max_B_femwell + +# # Inductance +# +# The inductance $L$ of a system given a current $I$ in a conductor can be calculated from +# +# $$ L = \frac{2W}{I^2} $$ +# +# with +# +# $$ W = \frac{1}{2} \int_\Omega B \cdot H d\Omega = \frac{1}{2} \int_\Omega \mu H \cdot H d\Omega $$ +# +# where $\mu$ is the permeability distribution, $B$ the magnetic field, $H = \mu^{-1} B$ the magnetization field, and $\Omega$ the domain. The integrand is only non-zero close to the conductors, where the field is concentrated. + +# We can compare the analytical solution to the FEM solution: diff --git a/femwell/magnetostatic.py b/femwell/magnetostatic.py new file mode 100644 index 00000000..0487a5bc --- /dev/null +++ b/femwell/magnetostatic.py @@ -0,0 +1,207 @@ +import matplotlib.pyplot as plt +import numpy as np +import shapely +from meshwell.model import Model +from meshwell.polysurface import PolySurface +from skfem import ( + Basis, + BilinearForm, + ElementDG, + ElementTriN1, + ElementTriP0, + ElementTriP1, + ElementVector, + FacetBasis, + Functional, + LinearForm, + asm, + condense, + enforce, + solve, +) +from skfem.helpers import curl, dot, grad, inner +from skfem.io import from_meshio +from skfem.io.meshio import from_meshio +from skfem.models.general import curluv +from tqdm import tqdm + +from femwell.coulomb import solve_coulomb +from femwell.mesh import mesh_from_OrderedDict +from femwell.visualization import plot_domains, plot_subdomain_boundaries + + +def solve_magnetostatic_2D(basis_mu_r, mu_r, mesh, current_densities): + """The 2D magnetostatic problem can be reduced to a scalar Coulomb equation with source.""" + + basis = basis_mu_r.with_element(ElementTriP1()) + + @BilinearForm + def dudv(A, v, w): + return w["mu_r"] * inner(grad(A), grad(v)) + + @LinearForm + def unit_load(v, _): + return v + + for ( + i, + (domain, current_density), + ) in enumerate(current_densities.items()): + basis_source = FacetBasis(mesh, basis.elem, facets=mesh.boundaries[domain], intorder=4) + basis_source_1d = basis_source.with_element(ElementTriP1()) + if i == 0: + J_rhs = basis_source_1d.zeros() + J_rhs += current_density * unit_load.assemble(basis_source_1d) + + A = dudv.assemble(basis, mu_r=basis_mu_r.interpolate(mu_r)) + A, J_rhs = enforce(A, J_rhs, D=mesh.boundary_nodes()) + + return basis, solve(A, J_rhs) + + +if __name__ == "__main__": + + # Define some parameters for the solenoid: + core_permeability = 1 + diameter = 2 + length = 1 + current = 1 + + def continuous_solenoid_cross_section_mesh( + diameter=diameter, + length=length, + dx_domain=10, + dy_domain=10, + ): + + solenoid_polygon = shapely.box(-diameter / 2, -length / 2, diameter / 2, length / 2) + left_air_polygon = shapely.box(-dx_domain / 2, -dy_domain / 2, -diameter / 2, dy_domain / 2) + right_air_polygon = shapely.box(diameter / 2, -dy_domain / 2, dx_domain / 2, dy_domain / 2) + center_air_polygon = shapely.box(-diameter / 2, -dy_domain / 2, diameter / 2, dy_domain / 2) + + model = Model() + + solenoid = PolySurface( + polygons=solenoid_polygon, + model=model, + physical_name="solenoid", + mesh_bool=True, + resolution={"resolution": 0.1, "DistMax": 2, "SizeMax": 0.2}, + mesh_order=1, + ) + left_air = PolySurface( + polygons=left_air_polygon, + model=model, + physical_name="left", + mesh_bool=True, + mesh_order=2, + ) + right_air = PolySurface( + polygons=right_air_polygon, + model=model, + physical_name="right", + mesh_bool=True, + mesh_order=2, + ) + center_air = PolySurface( + polygons=center_air_polygon, + model=model, + physical_name="center", + mesh_bool=True, + mesh_order=2, + ) + + return from_meshio( + model.mesh( + entities_list=[solenoid, left_air, center_air, right_air], + filename="mesh.msh", + default_characteristic_length=0.5, + verbosity=0, + ) + ) + + mesh = continuous_solenoid_cross_section_mesh( + diameter=diameter, + length=length, + dx_domain=10, + dy_domain=10, + ) + # mesh.draw().show() + + # plot_domains(mesh) + # plt.show() + + # plot_subdomain_boundaries(mesh) + + def magnetic_potential(mesh, current=1, permeability_mu_r=1): + + # Permeability defined constant-wise to handle sharp discontinuities + basis_mu_r = Basis(mesh, ElementTriP0()) + mu_r = basis_mu_r.ones() + mu_r[basis_mu_r.get_dofs(elements=("solenoid"))] = permeability_mu_r + + basis_A, A = solve_magnetostatic_2D( + basis_mu_r=basis_mu_r, + mu_r=mu_r, + mesh=mesh, + current_densities={ + "solenoid___left": current, + "solenoid___right": -current, + }, + ) + + return basis_A, A, basis_mu_r, mu_r + + basis_A, A, basis_mu_r, mu_r = magnetic_potential( + mesh, current=current, permeability_mu_r=core_permeability + ) + + fig, ax = plt.subplots() + for subdomain in ["solenoid"]: + basis_mu_r.mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True) + basis_A.plot(A, ax=ax, shading="gouraud", colorbar=True) + + bfield = asm( + LinearForm(curluv).partial(basis_A.interpolate(A)), + basis_A.with_element(ElementVector(ElementTriP1())), + ) + + print(np.min(basis_A.doflocs[0, :])) + print(np.max(basis_A.doflocs[0, :])) + print(np.min(basis_A.doflocs[1, :])) + print(np.max(basis_A.doflocs[1, :])) + + # Define the grid for interpolation + Nx, Ny = 50, 50 + xmin, xmax = mesh.p[0].min(), mesh.p[0].max() + ymin, ymax = mesh.p[1].min(), mesh.p[1].max() + grid_x, grid_y = np.meshgrid(np.linspace(xmin, xmax, Nx), np.linspace(ymin, ymax, Ny)) + + # Interpolate the magnetic field on the grid + from scipy.interpolate import griddata + + bfield_x = bfield.reshape((-1, 2))[:, 0] + bfield_y = bfield.reshape((-1, 2))[:, 1] + grid_bfield_x = griddata(mesh.p.T, bfield_x, (grid_x, grid_y), method="cubic") + grid_bfield_y = griddata(mesh.p.T, bfield_y, (grid_x, grid_y), method="cubic") + + # Plot the streamlines + ax.streamplot(grid_x, grid_y, grid_bfield_x, grid_bfield_y, color="black", linewidth=0.5) + + plt.show() + + # curl_A = curl(basis_A.interpolate(A)) + # print(A.shape) + # print(curl_A.shape) + # curl_A_x = basis_A.interpolate(curl_A[0,:]) + # curl_A_y = basis_A.interpolate(curl_A[1,:]) + # coordinates = basis_A.doflocs + + # velocity = asm( + # LinearForm(curluv).partial(ib.interpolate(psi)), + # ib.with_element(ElementVector(ElementTriP1())), + # ) + + # fig, ax = plt.subplots() + # ax.quiver(coordinates[0, :], coordinates[1, :], curl_A[0,:,0], curl_A[1,:,0]) + # plt.show() From 55da5f62a025ef1bd9fd1ebf4d6775dc4e612703 Mon Sep 17 00:00:00 2001 From: simbilod Date: Wed, 12 Jun 2024 23:42:14 -0700 Subject: [PATCH 3/4] update docs --- docs/_toc.yml | 2 + docs/math/maxwell.md | 35 -------- docs/math/quasistatic.md | 169 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 171 insertions(+), 35 deletions(-) create mode 100644 docs/math/quasistatic.md diff --git a/docs/_toc.yml b/docs/_toc.yml index 37029bf6..5bf6bf52 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -14,6 +14,7 @@ parts: - file: math/maxwell_quantities.md - file: math/waveguide_port.md - file: math/overlap_integrals.md + - file: math/quasistatic.md - file: math/thermal.md.md - file: math/dispersion.md - file: math/coupled_mode_theory.md @@ -45,6 +46,7 @@ parts: - file: photonics/examples/depletion_waveguide.py - file: photonics/examples/refinement.py - file: electronics/examples/capacitor.py + - file: electronics/examples/inductor.py - file: electronics/examples/coax_cable.py - file: electronics/examples/RF_CPW_transmission_line_tutorial/RF_CPW_transmission_line_tutorial.py - file: electronics/examples/coplanar_waveguide_vary_width.py diff --git a/docs/math/maxwell.md b/docs/math/maxwell.md index 348d5a3d..b5a7501d 100644 --- a/docs/math/maxwell.md +++ b/docs/math/maxwell.md @@ -301,41 +301,6 @@ where $R$ is the radius of curvature in $x$-direction. See discussion on choice of R in {cite}`Masi:10` -## Calculating static potentials - -As in the static case - -$$ - \nabla\times\vec{\mathcal{E}} - = - \mu \frac{\partial \vec{\mathcal{H}}}{\partial t} - = 0 -$$ - -$\mathcal{E}$ can be written as - -$$ - \vec{\mathcal{E}} = -\nabla \Phi -$$ (EdivPhi) - -using {eq}`maxwell`, for $\Phi$ can be found that - -$$ - -\nabla\cdot \left(\varepsilon \nabla \Phi\right) = \rho -$$ - -from which we can derive the weakform - -$$ - \left( - \varepsilon \nabla \Phi - , - \nabla v - \right) - = \rho v -$$ - -which is used to calculate the potential for a given structure. -Using {eq}`EdivPhi` the electric field can be calculated from the potential. ## 2D Periodic diff --git a/docs/math/quasistatic.md b/docs/math/quasistatic.md new file mode 100644 index 00000000..97858d06 --- /dev/null +++ b/docs/math/quasistatic.md @@ -0,0 +1,169 @@ +# Electromagnetic field simulations in the quasistatic case + +The below follows the treatment of \cite{Haus1989}, Ch 3. + +The quasistatic equations are derived by approximating one of the time derivatives in Maxwell's equations: +* $\frac{\partial \mathbf{B}\left(\mathbf{r},t\right)} {\partial t} \sim 0$: neglecting magnetic induction, leading to "electroquasistatic" (EQS) +* $\frac{\partial \mathbf{D}\left(\mathbf{r},t\right)} {\partial t} \sim 0$: neglecting electric displacement current, leading to "magnetoquasistatic" (MQS) + +These are justified if the condition + +$$ \frac{\mu \epsilon L^2}{\tau^2} \ll 1 \rightarrow \frac{L}{v} << \tau$$ + +i.e. if an electromagnetic wave can propagate a characteristic length of the system in a time much shorter than the timescales $\tau$ under consideration. For instance, with $v \sim c \approx 3 \times 10^{8}$ m/s, on the chip-scale with $L \sim 1 \times 10^{-3}$ m, the timescale considered much be much larger than $\sim 1 \times 10^{-12}$ s, in practical terms below GHz frequency regimes. + +Whether the electroquasistatic formulation or magnetoquasistatic formulation is used depends on which field dominates in the static limit. If the magnetic field vanishes, then the system is EQS (e.g. capacitor, non-looping conductor = resistor). Is the electric field vanishes, then the system is MQS (e.g. looping conductor = inductor). This is only easy to establish if the system is made up of "perfect conductors" and "perfect insulators", which at low frequencies regular metals and dielectrics are good approximations. + + +## Electroquasistatic case + +In the first case, we get irrotationality (curl-free) of the E-field + +$$ + \nabla\times\vec{\mathcal{E}} + = - \mu \frac{\partial \vec{\mathcal{H}}}{\partial t} + \sim 0 +$$ + +The lack of curl means that there is a scalar potential associated with $\mathcal{E}$: + +$$ + \vec{\mathcal{E}} = -\nabla \Phi +$$ (EdivPhi) + +This is the electric potential. When differences in the potential are referenced (to *e.g.* a "ground"), it is alternatively called voltage $V$. + +### In insulating systems + +This is explored more in-depth in \cite{Haus1989}, Ch. 4-6. + +using {eq}`maxwell`, for $\Phi$ can be found that + +$$ + -\nabla\cdot \left(\varepsilon \nabla \Phi\right) = \rho +$$ + +from which we can derive the weak form + +$$ + \left( + \varepsilon \nabla \Phi + , + \nabla v + \right) + = \rho v +$$ + +which is used to calculate the potential for a given structure. +Using {eq}`EdivPhi` the electric field can be calculated from the potential. + + +### In conducting systems + +This is explored more in-depth in \cite{Haus1989}, Ch. 7. + +Taking the divergence of Ampère–Maxwell's law, and substituting Gauss' law, we get + +$$ + \nabla \cdot \left( \nabla \times \vec{H} \right) = 0 =\nabla \cdot \left( \vec{J} + \frac{\partial \vec{D}}{\partial t} \right) + \rightarrow + \nabla \cdot \vec{J} = - \frac{\partial \rho}{\partial t} +$$ (continuity) + +This is also known as charge continuity. For DC currents (and without current sources), the assumption is that the net charge $\rho$ is unchanging, and hence $\nabla \cdot \vec{J} = 0$. + +Many materials exhibit a linear dependence between their supported current density and the applied voltage. This is known as Ohm's law: + +$$ \vec{J} = \sigma \vec{E} $$ + +Hence, in a purely conducting system, the electric field is the solution to + +$$ -\nabla \cdot \left( \sigma \nabla \Phi \right) = 0 $$ + +This has a similar weak form + +$$ + \left( + \sigma \nabla \Phi + , + \nabla v + \right) + = 0 +$$ + +### In combined systems + +When $\frac{\partial \rho}{\partial t}$ is nonzero, it can be included alongside $\vec{J}$ when taking the divergence. This leads to the displacement current $\frac{\partial \vec{D}}{\partial t}$ being added to the current: + +$$ + -\nabla \cdot \left( \sigma \nabla \Phi + \frac{\partial \left( \epsilon_0 \epsilon \nabla \Phi \right)}{\partial t} \right) = 0 +$$ + +or equivalently in Fourier space: + +$$ + -\nabla \cdot \left( \left( \sigma + i\omega \epsilon_o \epsilon \right) \nabla \Phi \right) = 0 +$$ + +This is exactly like the above with modified conduction. However, at DC ($\omega = 0$), in an inhomogeneous system, the current distribution is purely set by the relative conductivities. At low but nonzero frequencies, the permittivity acts as a frequency-dependent complex conductivity. + + +### Solving for H + +Once the potential (and hence E-field) has been calculated, the magnetic field can be calculated by using the remaining equations: + + + + +## Magnetoquasistatic case + +This is explored more in-depth in \cite{Haus1989}, Ch. 8-10, as well as \cite{Larson2013}, Ch. 13. + +In the second case, we get + +$$ + \nabla\times\vec{\mathcal{H}} + = \varepsilon \frac{\partial \vec{\mathcal{E}}}{\partial t} + \vec{J} + \sim \vec{J} +$$ (magnetic) + +combined with the no-monopole law $\nabla \cdot \vec{B} = 0$, we see that the magnetic field is solenoidal (curl, without divergence). + +The lack of divergence means that there is a vector potential associated with $\mathcal{B}$: + +$$ + \vec{\mathcal{B}} = \nabla \times \vec{A} +$$ (HcurlA) + +This is the magnetic vector potential. Previously in the electrostatic case we had a scalar potential, with arbitrarily-referenced potential $\Phi = 0$. Here we have a vector potential, with arbitrarily-referenced divergence $\nabla \cdot \vec{A} = 0$. This is the Coulomb gauge, chosen to simplify the below. + +Taking the curl of {eq}`magnetic`, substituting in the vector potential, and using $\nabla \cdot \vec{A} = 0$: + +$$ + \nabla \times \left( \mu^{-1} \nabla \times \vec{A} \right) = \vec{J} +$$ + +The first and easiest way to deal with this for piecewise constant permeability is to use the identity $\nabla \times \nabla \times \vec{V} = -\Delta \vec{V} + \nabla (\nabla \cdot \vec{V})$. In the Coulomb gauge, this yields a set of scalar equations for the components of $\vec A$. In 2D when $\vec{J} = J_z \hat{z}$, this is identical to the scalar Coulomb equation, solving for a scalar potential $\vec{A} = A_z \hat{z}$: + +$$ -\nabla \cdot \left( \mu^{-1} \nabla \cdot A_z \right) = J_z $$ + +In 3D, even though this form is in principle valid, the components remain coupled through boundary conditions. Hence, the full curl weak form with enforcement of 0 divergence through Lagrange multipliers is solved directly: + +$$\begin{aligned} + \left( \mu^{-1} \nabla \times \vec{A}, \nabla \times \vec{v} \right) - \left( v, \nabla \Lambda \right) - \left( \vec{A}, \nabla q \right) = \left( \vec{J}, \vec{v} \right) \\ +\end{aligned}$$ + +with $\vec{v} \in H_0(curl, \Omega)$ (curl-conforming), $q \in H_0^1(\Omega)$ (scalar), and $\left( \vec{A}, \Lambda \right) \in H_0(curl, \Omega) \times H_0^1(\Omega)$ + +which is used to calculate the potential for a given structure. +Using {eq}`HcurlA` the magnetic field can be calculated from the vector potential. + + + + +## Bibliography + +```{bibliography} +:style: unsrt +:filter: docname in docnames +``` From 5ab56c7a060820b14b8071ce673563ca78c79057 Mon Sep 17 00:00:00 2001 From: simbilod Date: Wed, 12 Jun 2024 23:48:56 -0700 Subject: [PATCH 4/4] add warning --- docs/electronics/examples/inductor.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/electronics/examples/inductor.py b/docs/electronics/examples/inductor.py index cb305fe4..8cc21312 100644 --- a/docs/electronics/examples/inductor.py +++ b/docs/electronics/examples/inductor.py @@ -16,6 +16,13 @@ # # We follow the theory outlined in the documentation to compute the magnetic field of various currents distributions. +#
+# The FEM results do not exactly yet match the analytical expectation: +# - There seems to be a sign error when calculating the magnetic field. +# - The potential values are slightly different. +# - The magnetic field values are very different. +#
+ # + import matplotlib.pyplot as plt import numpy as np