Heat Solver¶

This notebook demonstrates basic usage of Tidy3D's heat solver. We will consider three simple multimaterial setups, planar, cylindrical, and spherical, and compare results to analytical solutions.

Tidy3D's heat feature solves a steady-state heat equation: $$ - k \nabla^2 T = s $$ with temperature and flux continuity conditions on the interfaces between different materials: $$ T_1 = T_2 $$ $$ k_1 \frac{\partial T_1}{\partial n} = k_2 \frac{\partial T_2}{\partial n} $$ and subject to a set of boundary conditions on the simulation domain boundary consisting of the following types:

• Temperature (Dirichlet): $\displaystyle T = T_0$
• Heat flux (Neumann): $\displaystyle k \frac{\partial T}{\partial n} = f$
• Convection (Robin): $\displaystyle k \frac{\partial T}{\partial n} = c (T-T_0)$

Here, $k$ is the thermal conductivity, $s$ is the volumetric heat source, $T_0$ is the prescribed temperature value, $f$ is the prescribed heat flux, and $c$ is the convection heat transfer coefficient.

Package tidy3d contains all necessary API for defining and running heat simulation on our servers. Note that to use full visualization capabilities of heat solver data one needs package vtk. To install a compatible version of vtk one can use pip install tidy3d[vtk] command.

Defining Simulation Scene¶

First, we define basic parameter of the simulation setups. In all three cases, we will consider two material layers: one from coordinate r1 to r0, and the second one from r0 to r2. Thermal conductivities of the two layers are defined by parameters conductivity1 and conductivity2. Each layer contains a volumetric heat source with rate given by source1 and source2, correspondingly. Finally, the temperature at the bottom of the first layer is set to temperature1 while the temperature at the top of the second layer is temperature2. Parameters must be provided in units consistent with the spatial units used in Tidy3D, microns. In these tests we do not simulate some specific material, thus, vales of material properties are rather arbitrary.

To create materials for heat simulations, one must use heat_spec field of medium classes. Acceptable types for this field are SolidSpec and FluidSpec that define solid and fluid materials from the point of view of heat simulation. Note that the heat equation is only solved inside materials with SolidSpec, regions with FluidSpec or for which heat_spec is not provided at all are excluded in heat solver. Note that since currently Tidy3D supports only steady-state heat simulation parameter capacity will not have an effect on simulation results.

Let us first consider constructing structures for the planar case. We define three overlapping boxes, such that in the region from x=0 to x=1 we have a two layers consisting made of materials medium1 and medium2 on the background of background_medium.

Now we can define a Scene that specifies problem geometry for the planar case and visualize it. The purpose of class Scene is to allow visualization of simulation setups without providing physics-specific details such as sources, monitors, grid specifications, etc, and to provide an easy way to transfer simulation setups between different solvers.

Note that by default the region of plotting is determined automatically based on present in Scene structures. If all structures have infinite extent along a given axis (y in this particular case), the plotting limits are set to (-1, 1).

Custom plotting limits can be specified by using arguments vlim and hlim.

Analogously, we create three cylinders and three spheres in cylindrical and spherical cases, respectively. Note, for demonstration purposes we create one of the cylinders as an STL mesh using trimesh package.

Note that besides function .plot() that paints structures with colors not related to physical properties, there is also function .plot_heat_conductivity() that uses the heat conductivity of each structure to select a grayscale color for it. Note that regions where the heat equation is not solved (FluidSpec or heat_spec=None) are displayed using light blue hatching.

Defining Heat Simulation¶

Now we proceed to defining a specification for actual heat simulations. We begin with setting simulation domain sizes for each case. Since planar and cylindrical cases are essentially one- and two-dimensional setups, correspondingly, we set a small domain size along invariant directions.

Boundary Conditions¶

Boundary conditions for heat simulation are defined using class HeatBoundarySpec that has two required field condition and placement. The first one defines what boundary conditions should be imposed and accept types TemperatureBC, HeatFluxBC, or ConvectionBC. While the second one specifies where exactly it should be applied. Available options for that are:

• StructureBoundary: boundary of a structure. Only the portion of boundary not covered by subsequent structures is implied.
• StructureStructureInterface: interface between two structures. Specifically, the portion of boundary of the succeeding structure inside the preceding structure is implied.
• MediumMediumInterface: interface between two media.
• SimulationBoundary: boundary of heat simulation domain. Note that specific surfaces of simulation box can be chosen using field surfaces. For example: SimulationBoundary(surfaces=["z+", "z-"]). By default, all surfaces are selected.
• StructureSimulationBoundary: part of boundary of heat simulation domain covered by a structure. Again, only the portion of boundary not covered by subsequent structures is implied.

Note that to point at a particular structure or medium, one must set in their definitions field name and use it when creating a boundary condition placement.

For these simulations we set temperature boundary conditions on the top of upper layer and the bottom of the lower layer. We also explicitly set adiabatic boundary conditions on the simulation domain boundary. Note that any solution domain boundaries, which include simulation domain boundary and boundaries facing mediums with empty or FluidSpec type heat_spec, are assigned adiabatic boundary conditions if no other boundary conditions is specified. Thus, setting adiabatic boundary conditions on the simulation domain boundary is not necessary in these simulations and is only done for demonstration purposes.

Note also that, for example, in case of bc_bottom we could have equivalently set placement to be td.StructureBoundary(structure="core") or td.MediumMediumInterface(mediums=["solid1", "fluid"]).

Sources¶

We define volumetric heat source in each of the layers using HeatSource. Similarly to boundary conditions placement, heat sources are assigned to structures by using their names.

Monitors¶

Like in optic Tidy3D simulations, the results of calculations are stored using monitors. Currently, heat simulations support TemperatureMonitor that record spatial temperature distributions in chosen rectangular region. For convenience, even though the heat equation is solved on unstructured grids, TemperatureMonitor's can return simulation results on equivalent Cartesian grids in the form of SpatialDataArray. If the results on the original unstructured grid is desired, option unstructured=True can be used. In this case, the data will be presented in the form of TriangularGridDataset (for planar monitors) or TetrahedralGridDataset (for volumetric monitors). For more information about working with unstructured datasets please refer to the tutorial on unstructured data classes.

HeatChargeSimulation class¶

Combining purely geometric information about the simulation setup (Scene) and heat solver specific definitions we can create HeatChargeSimulation instance for each of the three cases. Structures and background medium of a Scene can either passed directly to a HeatChargeSimulation instance or a convenience method .from_scene() can be used. Note that in these examples we can make use of field symmetry to reduce the size and cost of simulations.

Plotting function of HeatChargeSimulation displays heat conductivity information of structures using gray scale colors, assigned boundary conditions using thick lines (orange for TemperatureBC, green for HeatFluxBC, and brown for ConvectionBC), and heat sources as colored dots. Specifying colorbar="conductivity" or colorbar="source" changes the colorbar accordingly.

Solving¶

Now we can submit the three heat simulation to our servers for solving. Heat simulations share the same web API as optic simulations. That is, we can submit simulations to servers in three ways.

The first way is to use function web.run(), which returns simulation results as a HeatChargeSimulationData instance.

Alternatively, one can use Job container.

Finally, class Batch facilitates a convenient submission of multiple simulations. Let us, for example, submit the spherical case at several different mesh resolutions.

For more details on Tidy3D's web API see the corresponding tutorial.

Data Visualization and Analysis¶

The solution data to the heat equation is stored in HeatChargeSimulationData object. Its field .data contains HeatMonitorData instances corresponding to each monitor from HeatChargeSimulation.monitor. Note that the heat equation is solved on an unstructured mesh, however, for a more convenient manipulation of data in the Tidy3D python client, the solution results can be transferred onto a Cartesian grid. If the results on the original unstructured grid is desired, option unstructured=True in monitor setup can be used. In that case, the data would be presented in the form of TriangularGridDataset (for planar monitors) or TetrahedralGridDataset (for volumetric monitors). For more information about working with unstructured datasets please refer to the tutorial on unstructured data classes.

Temperature distribution corresponding to a given monitor can be visualized using function HeatChargeSimulationData.plot_field(). Note that numerical solutions are slightly extrapolated outside of solution domains, such that interpolation can be accurately performed anywhere in the vicinity of the solution domain. Additionally, data in monitor data objects are automatically expanded with the respect to the symmetry of the original simulation.

Now let us compare the numerical results to analytical solutions. The general solution to the quasi one-dimensional steady-state heat equation in 1, 2 and 3 dimensions can be written as $$ T_a(r) = - \frac1{2N} \frac{s}{k} r^2 + A u_N(r) + B,$$ where $N$ is the problem's dimensionality, $s$ is the source term, $k$ is the heat conductivity, $A$ and $B$ are coefficients defined by boundary conditions, and $u_N(r)$ is the Green function, which is given by $$ u_1(r) = r, \quad u_2(r) = \log(r), \quad u_3(r) = - \frac{1}{r}, $$ in one, two, and three dimensions, respectively.

The analytical solutions $T(r)$ for the considered two-layer setups can be decomposed into two parts $T_1(r)$ and $T_2(r)$ such that $$ T(r) = T_1(r),\, r < r_0$$ $$ T(r) = T_2(r),\, r > r_0$$ and each of which are given by the form $T_a(r)$ but with different coefficients $A_1$, $B_1$ and $A_2$, $B_2$. These are defined from boundary conditions: $$T_1(r_1) = t_1,$$ $$T_2(r_2) = t_2,$$ $$T_1(r_0) = T_2(r_0),$$ $$k_1 \frac{d}{dr} T_1(r_0) = k_2\frac{d}{dr} T_2(r_0).$$ Thus, coefficients $A_1$, $A_2$, $B_1$, and $B_2$ can be found by solving linear system $$ \begin{pmatrix} u_N(r_1)& 1& 0& 0 \\ 0& 0& u_N(r_2)& 1 \\ u_N(r_0)& 1& -u(r_0)& -1 \\ k_1 \frac{du_N}{dr}(r_0)& 0& -k_2 * \frac{du_N}{dr}(r_0)& 0 \end{pmatrix} \begin{pmatrix} A_1 \\ B_1 \\ A_2 \\ B_2 \end{pmatrix} = \begin{pmatrix} t_1 + \frac1{2N} \frac{s_1}{k_1} r_1^2 \\ t_2 + \frac1{2N} \frac{s_2}{k_2} r_2^2 \\ \frac1{2N} r_0^2 \left(\frac{s_1}{k_1} - \frac{s_2}{k_2}\right) \\ \frac1{N} r_0 (s_1 - s_2), \end{pmatrix} $$

Let us define a helper function that finds coefficients $A$ and $B$ from given boundary conditions.

And another helper function that samples analytical solution at given coordinate points.

Now we can easily find analytical solution for each of the cases and compare them to the numerical ones. Let us first perform this comparison for all three cases at resolution dl = 0.02.

As one can see the solver accurately reproduces analytical solutions.

Now, let us look at the convergence of numerical solutions in the spherical case. Plotting the dependence of the solution error (in the $L_\infty$ norm) on the grid resolution we confirm the second order of convergence of the solver.