Radiative cooling glass coating#
Air conditioning accounts for about 10% of global electricity use, and demand is projected to rise significantly. New approaches like radiative cooling can reduce this energy demand and associated greenhouse gas emissions. Radiative cooling uses materials applied to buildings that reflect sunlight while emitting heat through the atmospheric transparency window (8-13 μm) to outer space (~3K). This can passively cool buildings below ambient temperature. An ideal radiative cooling material would combine high solar reflection and selective infrared emission for cooling, easy manufacturability at scale, mechanical robustness, and long-term stability under outdoor environmental exposure. A recent paper introduces a new stable, scalable, and low-cost radiative cooling coating using a microporous glass-Al\(_2\)O\(_3\) particle composite to address these needs.
This notebook simulates the emissivity of the 100 μm thick coating and the result shows that the emissivity (absorptance) of the coating is near unity in the atmospheric transparency window. The design is based on Xinpeng Zhao et al., A solution-processed radiative cooling glass. Science 382, 684-691(2023).DOI: 10.1126/science.adi2224.
This notebook is contributed by Dr. Yurui Qu.
[1]:
import matplotlib.pylab as plt
import numpy as np
import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins.dispersion import AdvancedFastFitterParam, FastDispersionFitter
Simulation Setup#
The 100 μm thick coating consists of randomly placed Al\(_2\)O\(_3\) spheres with 250 nm radius and SiO\(_2\) spheres with 4 μm radius. The Al\(_2\)O\(_3\) spheres and the and SiO\(_2\) spheres take up 20% and 30% of the total volume of the coating layer, respectively. We truncate the coating to 20 μm in the in-plane directions and will apply the periodic boundary condition.
[2]:
# Define Parameters
# radius and location of the sphere
radius_Al2O3 = 0.25
radius_SiO2 = 4 # exp is 6um
box_size_xy = 20
box_size_z = 100
vol_Al2O3 = 4 / 3 * np.pi * np.power(radius_Al2O3, 3)
vol_SiO2 = 4 / 3 * np.pi * np.power(radius_SiO2, 3)
vol_box = box_size_xy * box_size_xy * box_size_z
num_Al2O3 = int(np.floor(0.2 * vol_box / vol_Al2O3)) # 20% of volume is Al2O3
num_SiO2 = int(np.floor(0.3 * vol_box / vol_SiO2)) # 30% of volume is SiO2
print("num_Al2O3:", num_Al2O3)
print("num_SiO2:", num_SiO2)
num_Al2O3: 122230
num_SiO2: 44
The refractive indices of Al\(_2\)O\(_3\) and SiO\(_2\) will be read from .csv files and fit with Tidy3D’s FastDispersionFitter. Due to the phonon resonances in the mid-infrared, both materials require a relatively large number of poles to yield a good fit.
[3]:
# permittivity of Al2O3
mat_Al2O3 = "misc/mat_Al2O3.csv"
advanced_param = AdvancedFastFitterParam(weights=(1, 1))
fitter = FastDispersionFitter.from_file(mat_Al2O3, skiprows=1, delimiter=",")
medium_Al2O3, rms_error = fitter.fit(
max_num_poles=6, advanced_param=advanced_param, tolerance_rms=2e-2
)
fitter.plot(medium_Al2O3)
plt.show()
14:25:34 CEST WARNING: Unable to fit with weighted RMS error under 'tolerance_rms' of 0.02
[4]:
mat_SiO2 = "misc/mat_SiO2.csv"
fitter = FastDispersionFitter.from_file(mat_SiO2, skiprows=1, delimiter=",")
medium_SiO2, rms_error = fitter.fit(
max_num_poles=8, advanced_param=advanced_param, tolerance_rms=2e-2
)
fitter.plot(medium_SiO2)
plt.show()
14:30:52 CEST WARNING: Unable to fit with weighted RMS error under 'tolerance_rms' of 0.02
Define some basic simulation parameters.
[5]:
# free space central wavelength
wl_start = 4 # wavelength
wl_end = 20 # wavelength
freq_start = td.C_0 / wl_end
freq_end = td.C_0 / wl_start
freqs = np.linspace(freq_start, freq_end, 100) # frequency range of the simulation
freq0 = (freq_start + freq_end) / 2 # central frequency
freqw = freq_end - freq_start # width of the frequency range
# distance between the surface of the sphere and the start of the PML layers along each cartesian direction
buffer_PML = 2 * wl_end
buffer_source = 1 * wl_end
# set the domain size in x, y, and z
domain_size_xy = box_size_xy
domain_size_z = buffer_PML + box_size_z + buffer_PML
# construct simulation size array
sim_size = (domain_size_xy, domain_size_xy, domain_size_z)
Next, we generate the randomly placed spheres.
[6]:
# Create random structures
SiO2_geometry = []
Al2O3_geometry = []
geometry = []
for i in range(num_SiO2):
position_xy = (box_size_xy - 2 * radius_SiO2) * (np.random.rand(2) - 0.5)
position_z = (box_size_z - 2 * radius_SiO2) * (np.random.rand(1) - 0.5)
position = [position_xy[0], position_xy[1], position_z]
sphere = td.Sphere(center=position, radius=radius_SiO2)
SiO2_geometry.append(sphere)
geometry.append(
td.Structure(geometry=td.GeometryGroup(geometries=SiO2_geometry), medium=medium_SiO2)
)
for i in range(num_Al2O3):
position_xy = (box_size_xy - 2 * radius_Al2O3) * (np.random.rand(2) - 0.5)
position_z = (box_size_z - 2 * radius_Al2O3) * (np.random.rand(1) - 0.5)
position = [position_xy[0], position_xy[1], position_z]
sphere = td.Sphere(center=position, radius=radius_Al2O3)
Al2O3_geometry.append(sphere)
geometry.append(
td.Structure(geometry=td.GeometryGroup(geometries=Al2O3_geometry), medium=medium_Al2O3)
)
geometry = tuple(geometry)
A PlaneWave is defined above the coating layer and two FluxMonitors are defined on both sides of the coating layer to measure transmission and reflection. Lastly, a FieldMonitor is defined to help visualize the field distribution within the coating layer.
[7]:
# add a plane wave source
plane_wave = td.PlaneWave(
source_time=td.GaussianPulse(freq0=freq0, fwidth=0.5 * freqw),
size=(td.inf, td.inf, 0),
center=(0, 0, box_size_z / 2 + buffer_source),
direction="-",
pol_angle=0,
)
# add a flux monitor to detect transmission
monitor_t = td.FluxMonitor(
center=[0, 0, -box_size_z / 2 - (buffer_source + buffer_PML) / 2],
size=[td.inf, td.inf, 0],
freqs=freqs,
name="T",
)
# add a flux monitor to detect reflection
monitor_r = td.FluxMonitor(
center=[0, 0, box_size_z / 2 + (buffer_source + buffer_PML) / 2],
size=[td.inf, td.inf, 0],
freqs=freqs,
name="R",
)
# add a field monitor to see the field profile at the absorption peak frequency
monitor_field = td.FieldMonitor(
center=[0, 0, 0], size=[td.inf, 0, td.inf], freqs=[freq0], name="field"
)
Define a Tidy3D Simulation.
[8]:
run_time = 2e-11 # simulation run time
# set up simulation
sim = td.Simulation(
size=sim_size,
grid_spec=td.GridSpec.uniform(dl=wl_start / 20),
structures=geometry,
sources=[plane_wave],
monitors=[monitor_t, monitor_r, monitor_field],
run_time=run_time,
boundary_spec=td.BoundarySpec(
x=td.Boundary.periodic(), y=td.Boundary.periodic(), z=td.Boundary.pml()
),
)
14:31:13 CEST WARNING: Structure at 'structures[1]' has bounds that extend exactly to simulation edges. This can cause unexpected behavior. If intending to extend the structure to infinity along one dimension, use td.inf as a size variable instead to make this explicit.
Before running the simulation, we can inspect a few things to ensure the simulation is set up correctly. First, plot the source time and frequency spectrum.
[9]:
# Visualize source
plane_wave.source_time.plot(np.linspace(0, run_time / 10, 1001))
plt.show()
plane_wave.source_time.plot_spectrum(times=np.linspace(0, run_time / 10, 2000), val="abs")
plt.show()
Then, plot the cross sectional view of the simulation.
[10]:
cfig, ax = plt.subplots(1, 2, figsize=(7, 3))
sim.plot(y=0, ax=ax[0])
sim.plot(z=0.1, freq=freq0, ax=ax[1])
plt.show()
Running the Simulation#
Once we confirm that the simulation is set up correctly, we can upload the simulation and calculate the maximum FlexCredit cost. This step prevents us from submitting large simulations by mistake.
[11]:
task_id = web.upload(sim, task_name="Simulation")
estimated_cost = web.estimate_cost(task_id)
print(f"The estimated maximum cost is {estimated_cost:.3f} Flex Credits.")
14:32:37 CEST Created task 'Simulation' with task_id 'fdve-21f61fb5-d8bc-46df-a6ee-fd3d373c7ed0' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-21f61fb5-d8 bc-46df-a6ee-fd3d373c7ed0'.
Task folder: 'default'.
14:33:52 CEST Maximum FlexCredit cost: 0.725. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
14:34:35 CEST Maximum FlexCredit cost: 0.725. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
The estimated maximum cost is 0.725 Flex Credits.
The cost seems reasonably so we start the task and monitor its status. After the simulation is complete, we can print out the real cost.
[ ]:
web.start(task_id)
web.monitor(task_id, verbose=True)
import time
time.sleep(20)
print("Billed flex unit cost: ", web.real_cost(task_id))
sim_data = web.load(task_id, path="data/cooling.hdf5")
# Show the output of the log file
print(sim_data.log)
14:34:36 CEST status = queued
To cancel the simulation, use 'web.abort(task_id)' or 'web.delete(task_id)' or abort/delete the task in the web UI. Terminating the Python script will not stop the job running on the cloud.
14:34:55 CEST status = preprocess
14:35:04 CEST starting up solver
running solver
14:35:45 CEST early shutoff detected at 36%, exiting.
14:35:46 CEST status = postprocess
14:36:07 CEST status = success
14:36:09 CEST View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-21f61fb5-d8 bc-46df-a6ee-fd3d373c7ed0'.
14:36:29 CEST Billed flex credit cost: 0.261.
Note: the task cost pro-rated due to early shutoff was below the minimum threshold, due to fast shutoff. Decreasing the simulation 'run_time' should decrease the estimated, and correspondingly the billed cost of such tasks.
Billed flex unit cost: 0.2610674426389001
14:36:41 CEST loading simulation from data/cooling.hdf5
14:37:08 CEST WARNING: Structure at 'structures[1]' has bounds that extend exactly to simulation edges. This can cause unexpected behavior. If intending to extend the structure to infinity along one dimension, use td.inf as a size variable instead to make this explicit.
Result Visualization#
After the simulation is complete, we will plot the transmission, reflection, and absorption spectra. The result shows that the absorption is nearly 100% in the atmospheric transparency window of 8-13 μm. Furthermore, based on Kirchhoff’s law of thermal radiation, the emissivity of a material is equal to its absorptance, we know that the infrared emissivity of the coating layer is close to unity.
Lastly, we can plot the field distribution within the coating layer. The plot shows that the electromagnetic energy is strongly absorbed as it propagates into the layer.
[ ]:
# Result Visualization
R = sim_data["R"].flux
T = -sim_data["T"].flux
A = 1 - R - T
plt.plot(td.C_0 / freqs, R, td.C_0 / freqs, T, td.C_0 / freqs, A)
# Save the absorption spectrum as as a .txt file
np.savetxt("data/Abs_4-20um.txt", (np.transpose((td.C_0 / freqs, A))))
plt.xlabel("Wavelength (μm)")
plt.ylim(0, 1)
plt.legend(("R", "T", "A"))
plt.show()
ax = sim_data.plot_field(field_monitor_name="field", field_name="E", val="abs^2")
ax.set_aspect("auto")
plt.show()
[ ]: