heater notebook

docs/notebooks/004_heater.py

@@ -0,0 +1,187 @@
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.15.2
+#   kernelspec:
+#     display_name: Python 3
+#     name: python3
+# ---
+
+# # TiN TOPS heater
+
+# + tags=["hide-input"]
+from collections import OrderedDict
+
+import matplotlib.pyplot as plt
+import numpy as np
+from femwell.maxwell.waveguide import compute_modes
+from femwell.mesh import mesh_from_OrderedDict
+from femwell.thermal import solve_thermal
+from shapely.geometry import LineString, Polygon
+from skfem import Basis, ElementTriP0
+from skfem.io import from_meshio
+from tqdm import tqdm
+
+# -
+
+# Simulating the TiN TOPS heater in {cite}`Jacques2019`.
+# First we set up the mesh:
+
+# + tags=["remove-stderr"]
+
+w_sim = 40  # 8 * 2
+h_clad = 2.8
+h_box = 2
+w_core = 40
+h_core = 0.22
+h_heater = 0.14
+w_heater = 2
+offset_heater = 2 + (h_core + h_heater) / 2
+h_silicon = 0.5
+
+polygons = OrderedDict(
+    bottom=LineString(
+        [
+            (-w_sim / 2, -h_core / 2 - h_box - h_silicon),
+            (w_sim / 2, -h_core / 2 - h_box - h_silicon),
+        ]
+    ),
+    core=Polygon(
+        [
+            (-w_core / 2, -h_core / 2),
+            (-w_core / 2, h_core / 2),
+            (w_core / 2, h_core / 2),
+            (w_core / 2, -h_core / 2),
+        ]
+    ),
+    heater=Polygon(
+        [
+            (-w_heater / 2, -h_heater / 2 + offset_heater),
+            (-w_heater / 2, h_heater / 2 + offset_heater),
+            (w_heater / 2, h_heater / 2 + offset_heater),
+            (w_heater / 2, -h_heater / 2 + offset_heater),
+        ]
+    ),
+    clad=Polygon(
+        [
+            (-w_sim / 2, -h_core / 2),
+            (-w_sim / 2, -h_core / 2 + h_clad),
+            (w_sim / 2, -h_core / 2 + h_clad),
+            (w_sim / 2, -h_core / 2),
+        ]
+    ),
+    box=Polygon(
+        [
+            (-w_sim / 2, -h_core / 2),
+            (-w_sim / 2, -h_core / 2 - h_box),
+            (w_sim / 2, -h_core / 2 - h_box),
+            (w_sim / 2, -h_core / 2),
+        ]
+    ),
+    wafer=Polygon(
+        [
+            (-w_sim / 2, -h_core / 2 - h_box - h_silicon),
+            (-w_sim / 2, -h_core / 2 - h_box),
+            (w_sim / 2, -h_core / 2 - h_box),
+            (w_sim / 2, -h_core / 2 - h_box - h_silicon),
+        ]
+    ),
+)
+
+resolutions = dict(
+    # core={"resolution": 0.04, "distance": 1},
+    clad={"resolution": 0.6, "distance": 1},
+    box={"resolution": 0.6, "distance": 1},
+    heater={"resolution": 0.1, "distance": 1},
+)
+
+mesh = from_meshio(
+    mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6)
+)
+mesh.draw().show()
+# -
+
+# And then we solve it!
+
+# + tags=["remove-stderr"]
+currents = np.linspace(0.0, 7.4e-3, 10)
+current_densities = currents / polygons["heater"].area
+neffs = []
+
+for current_density in tqdm(current_densities):
+    basis0 = Basis(mesh, ElementTriP0(), intorder=4)
+    thermal_conductivity_p0 = basis0.zeros()
+    for domain, value in {
+        "core": 90,
+        "box": 1.38,
+        "clad": 1.38,
+        "heater": 28,
+        "wafer": 148,
+    }.items():
+        thermal_conductivity_p0[basis0.get_dofs(elements=domain)] = value
+    thermal_conductivity_p0 *= 1e-12  # 1e-12 -> conversion from 1/m^2 -> 1/um^2
+
+    basis, temperature = solve_thermal(
+        basis0,
+        thermal_conductivity_p0,
+        specific_conductivity={"heater": 2.3e6},
+        current_densities={"heater": current_density},
+        fixed_boundaries={"bottom": 0},
+    )
+
+    if current_density == current_densities[-1]:
+        basis.plot(temperature, shading="gouraud", colorbar=True)
+        plt.show()
+
+    temperature0 = basis0.project(basis.interpolate(temperature))
+    epsilon = basis0.zeros() + (1.444 + 1.00e-5 * temperature0) ** 2
+    epsilon[basis0.get_dofs(elements="core")] = (
+        3.4777 + 1.86e-4 * temperature0[basis0.get_dofs(elements="core")]
+    ) ** 2
+    # basis0.plot(epsilon, colorbar=True).show()
+
+    modes = compute_modes(basis0, epsilon, wavelength=1.55, num_modes=1)
+
+    if current_density == current_densities[-1]:
+        modes[0].show(modes[0].E.real)
+
+    neffs.append(np.real(modes[0].n_eff))
+
+print(f"Phase shift: {2 * np.pi / 1.55 * (neffs[-1] - neffs[0]) * 320}")
+plt.xlabel("Current / mA")
+plt.ylabel("Effective refractive index $n_{eff}$")
+plt.plot(currents * 1e3, neffs)
+plt.show()
+
+
+# +
+dofs = basis.get_dofs(elements="core").flatten()
+
+plt.scatter(
+    basis.doflocs[:, dofs.flatten()][0, :],
+    temperature[dofs] / np.max(temperature[dofs]),
+    marker=".",
+)
+plt.axvline(x=0, color="b", linestyle="--")
+plt.axvline(x=2, color="g", linestyle="--")
+plt.axvline(x=5, color="r", linestyle="--")
+plt.axvline(x=8, color="c", linestyle="--")
+plt.axvline(x=11, color="m", linestyle="--")
+plt.axvline(x=20, color="y", linestyle="--")
+
+plt.title("Diffusion length at core layer")
+plt.xlabel("x (um)")
+plt.ylabel("Temperature (relative to underneath heater)")
+
+# -
+
+# ## Bibliography
+#
+# ```{bibliography}
+# :style: unsrt
+# :filter: docname in docnames
+# ```