From a058736d3ec7cb01ab804a69897248ce7d01d272 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Damiano=20Franz=C3=B2?= Date: Thu, 24 Jul 2025 17:53:19 +0200 Subject: [PATCH 1/9] FEAT: Add natural convection BC --- tidy3d/__init__.py | 3 + tidy3d/components/tcad/boundary/heat.py | 101 +++++++++++++++++++++++- tidy3d/components/tcad/types.py | 2 +- 3 files changed, 104 insertions(+), 2 deletions(-) diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index a3993e3863..f7bf678566 100644 --- a/tidy3d/__init__.py +++ b/tidy3d/__init__.py @@ -68,6 +68,7 @@ CaugheyThomasMobility, ConstantMobilityModel, ConvectionBC, + NaturalConvectionVerticalSpec, CurrentBC, HeatFluxBC, HeatFromElectricSource, @@ -405,6 +406,7 @@ def set_logging_level(level: str) -> None: Transformed.update_forward_refs() ClipOperation.update_forward_refs() GeometryGroup.update_forward_refs() +ConvectionBC.update_forward_refs() __all__ = [ "C_0", @@ -713,4 +715,5 @@ def set_logging_level(level: str) -> None: "set_logging_console", "set_logging_file", "wavelengths", + "NaturalConvectionVerticalSpec", ] diff --git a/tidy3d/components/tcad/boundary/heat.py b/tidy3d/components/tcad/boundary/heat.py index 9671dee6a9..e9fc097b59 100644 --- a/tidy3d/components/tcad/boundary/heat.py +++ b/tidy3d/components/tcad/boundary/heat.py @@ -4,6 +4,7 @@ import pydantic.v1 as pd +from typing import Union from tidy3d.components.tcad.boundary.abstract import HeatChargeBC from tidy3d.constants import HEAT_FLUX, HEAT_TRANSFER_COEFF, KELVIN @@ -39,6 +40,100 @@ class HeatFluxBC(HeatChargeBC): units=HEAT_FLUX, ) +#TODO Should I adapt to Tidy3D units (e.g. micrometers), or keep the original ones? +THERMAL_CONDUCTIVITY_UNITS = "W/(m*K)" +DYNAMIC_VISCOSITY_UNITS = "Pa*s" +SPECIFIC_HEAT_UNITS = "J/(kg*K)" +DENSITY_UNITS = "kg/m**3" +THERMAL_EXPANSIVITY_UNITS = "1/K" +LENGTH_UNITS = "m" +ACCELERATION_UNITS = "m/s**2" + +class NaturalConvectionVerticalSpec(HeatChargeBC): + """ + Specification for natural convection from a vertical plate. + + This class calculates the heat transfer coefficient (h) based on fluid + properties and an expected temperature difference, then provides these + values as 'base' and 'exponent' for a generalized heat flux equation + q = base * (T_surf - T_fluid)^exponent. + """ + # --- Input Parameters --- + fluid_k: pd.NonNegativeFloat = pd.Field( + title="Fluid Thermal Conductivity", + description="Thermal conductivity (k) of the fluid.", + units=THERMAL_CONDUCTIVITY_UNITS, + ) + fluid_mu: pd.NonNegativeFloat = pd.Field( + title="Fluid Dynamic Viscosity", + description="Dynamic viscosity (μ) of the fluid.", + units=DYNAMIC_VISCOSITY_UNITS, + ) + fluid_Cp: pd.NonNegativeFloat = pd.Field( + title="Fluid Specific Heat", + description="Specific heat capacity (Cp) of the fluid at constant pressure.", + units=SPECIFIC_HEAT_UNITS, + ) + fluid_rho: pd.NonNegativeFloat = pd.Field( + title="Fluid Density", + description="Density (ρ) of the fluid.", + units=DENSITY_UNITS, + ) + fluid_beta: pd.NonNegativeFloat = pd.Field( + title="Fluid Thermal Expansivity", + description="Thermal expansion coefficient (β) of the fluid.", + units=THERMAL_EXPANSIVITY_UNITS, + ) + plate_L: pd.NonNegativeFloat = pd.Field( + title="Plate Characteristic Length", + description="Characteristic length (L), defined as the height of the vertical plate.", + units=LENGTH_UNITS, + ) + + gravity: pd.NonNegativeFloat = pd.Field( + default=9.80665, + title="Gravitational Acceleration", + description="Gravitational acceleration (g).", + units=ACCELERATION_UNITS, + ) + + def _compute_parameters(self): + + # Calculate the Rayleigh Number (Ra_L) + rayleigh_numerator_notemp = ( + self.gravity + * self.fluid_beta + * self.fluid_rho ** 2 + * self.fluid_Cp + * self.plate_L ** 3 + ) + rayleigh_denominator = self.fluid_mu * self.fluid_k + Ra_L = rayleigh_numerator_notemp / rayleigh_denominator + + # Calculate the denominator term from the Nusselt number correlation + # This term is related to the Prandtl Number (Pr = mu * Cp / k) + pr_term_inner_num = 0.492 * self.fluid_k + pr_term_inner_den = self.fluid_mu * self.fluid_Cp + + pr_term_inner = pr_term_inner_num / pr_term_inner_den + pr_denominator = (1 + pr_term_inner ** (9 / 16)) ** (4 / 9) + + # Select formula based on flow regime (determined by Ra_L) + if Ra_L <= 1e9: + # Laminar Flow + h_factor_linear = 0.68 + h_factor_non_linear = (0.670 * Ra_L ** (1 / 6)) / pr_denominator + elif 1e9 < Ra_L <= 1e13: + # Turbulent Flow + h_factor_linear = 0.825 + h_factor_non_linear = (0.387 * Ra_L ** (1 / 6)) / pr_denominator + else: + raise ValueError(f"Ra_l={Ra_L} should be smaller than 1e13 for NaturalConvectionVerticalSpec") + + h = (self.fluid_k / self.plate_L) * h_factor_linear + h_nonlinear = (self.fluid_k / self.plate_L) * h_factor_non_linear + exponent = 1 + 1/6 + return h, h_nonlinear, exponent class ConvectionBC(HeatChargeBC): """Convective thermal boundary conditions. @@ -55,8 +150,12 @@ class ConvectionBC(HeatChargeBC): units=KELVIN, ) - transfer_coeff: pd.NonNegativeFloat = pd.Field( + transfer_coeff: Union[pd.NonNegativeFloat, NaturalConvectionVerticalSpec] = pd.Field( title="Heat Transfer Coefficient", description=f"Heat flux value in units of {HEAT_TRANSFER_COEFF}.", units=HEAT_TRANSFER_COEFF, ) + +#TODO Maybe I should find a better place for these lines. +NaturalConvectionVerticalSpec.update_forward_refs() +ConvectionBC.update_forward_refs() diff --git a/tidy3d/components/tcad/types.py b/tidy3d/components/tcad/types.py index dec3f44f62..79fd754f55 100644 --- a/tidy3d/components/tcad/types.py +++ b/tidy3d/components/tcad/types.py @@ -4,7 +4,7 @@ from tidy3d.components.tcad.bandgap import SlotboomBandGapNarrowing from tidy3d.components.tcad.boundary.charge import CurrentBC, InsulatingBC, VoltageBC -from tidy3d.components.tcad.boundary.heat import ConvectionBC, HeatFluxBC, TemperatureBC +from tidy3d.components.tcad.boundary.heat import ConvectionBC, HeatFluxBC, TemperatureBC, NaturalConvectionVerticalSpec from tidy3d.components.tcad.generation_recombination import ( AugerRecombination, RadiativeRecombination, From bbe0d31e5189236d2f0e6a7629199a3d9c482858 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Damiano=20Franz=C3=B2?= Date: Thu, 24 Jul 2025 17:59:57 +0200 Subject: [PATCH 2/9] Fixing formatting and minor error --- tidy3d/components/tcad/boundary/heat.py | 24 +++++++++++++----------- tidy3d/components/tcad/types.py | 7 ++++++- 2 files changed, 19 insertions(+), 12 deletions(-) diff --git a/tidy3d/components/tcad/boundary/heat.py b/tidy3d/components/tcad/boundary/heat.py index e9fc097b59..22e3567691 100644 --- a/tidy3d/components/tcad/boundary/heat.py +++ b/tidy3d/components/tcad/boundary/heat.py @@ -40,7 +40,8 @@ class HeatFluxBC(HeatChargeBC): units=HEAT_FLUX, ) -#TODO Should I adapt to Tidy3D units (e.g. micrometers), or keep the original ones? + +# TODO Should I adapt to Tidy3D units (e.g. micrometers), or keep the original ones? THERMAL_CONDUCTIVITY_UNITS = "W/(m*K)" DYNAMIC_VISCOSITY_UNITS = "Pa*s" SPECIFIC_HEAT_UNITS = "J/(kg*K)" @@ -49,6 +50,7 @@ class HeatFluxBC(HeatChargeBC): LENGTH_UNITS = "m" ACCELERATION_UNITS = "m/s**2" + class NaturalConvectionVerticalSpec(HeatChargeBC): """ Specification for natural convection from a vertical plate. @@ -58,6 +60,7 @@ class NaturalConvectionVerticalSpec(HeatChargeBC): values as 'base' and 'exponent' for a generalized heat flux equation q = base * (T_surf - T_fluid)^exponent. """ + # --- Input Parameters --- fluid_k: pd.NonNegativeFloat = pd.Field( title="Fluid Thermal Conductivity", @@ -98,14 +101,9 @@ class NaturalConvectionVerticalSpec(HeatChargeBC): ) def _compute_parameters(self): - # Calculate the Rayleigh Number (Ra_L) rayleigh_numerator_notemp = ( - self.gravity - * self.fluid_beta - * self.fluid_rho ** 2 - * self.fluid_Cp - * self.plate_L ** 3 + self.gravity * self.fluid_beta * self.fluid_rho**2 * self.fluid_Cp * self.plate_L**3 ) rayleigh_denominator = self.fluid_mu * self.fluid_k Ra_L = rayleigh_numerator_notemp / rayleigh_denominator @@ -122,19 +120,22 @@ def _compute_parameters(self): if Ra_L <= 1e9: # Laminar Flow h_factor_linear = 0.68 - h_factor_non_linear = (0.670 * Ra_L ** (1 / 6)) / pr_denominator + h_factor_non_linear = (0.670 * Ra_L ** (1 / 6)) / pr_denominator elif 1e9 < Ra_L <= 1e13: # Turbulent Flow h_factor_linear = 0.825 h_factor_non_linear = (0.387 * Ra_L ** (1 / 6)) / pr_denominator else: - raise ValueError(f"Ra_l={Ra_L} should be smaller than 1e13 for NaturalConvectionVerticalSpec") + raise ValueError( + f"Ra_L={Ra_L} should be smaller than 1e13 for NaturalConvectionVerticalSpec" + ) h = (self.fluid_k / self.plate_L) * h_factor_linear h_nonlinear = (self.fluid_k / self.plate_L) * h_factor_non_linear - exponent = 1 + 1/6 + exponent = 1 + 1 / 6 return h, h_nonlinear, exponent + class ConvectionBC(HeatChargeBC): """Convective thermal boundary conditions. @@ -156,6 +157,7 @@ class ConvectionBC(HeatChargeBC): units=HEAT_TRANSFER_COEFF, ) -#TODO Maybe I should find a better place for these lines. + +# TODO Maybe I should find a better place for these lines. NaturalConvectionVerticalSpec.update_forward_refs() ConvectionBC.update_forward_refs() diff --git a/tidy3d/components/tcad/types.py b/tidy3d/components/tcad/types.py index 79fd754f55..38cb64c74b 100644 --- a/tidy3d/components/tcad/types.py +++ b/tidy3d/components/tcad/types.py @@ -4,7 +4,12 @@ from tidy3d.components.tcad.bandgap import SlotboomBandGapNarrowing from tidy3d.components.tcad.boundary.charge import CurrentBC, InsulatingBC, VoltageBC -from tidy3d.components.tcad.boundary.heat import ConvectionBC, HeatFluxBC, TemperatureBC, NaturalConvectionVerticalSpec +from tidy3d.components.tcad.boundary.heat import ( + ConvectionBC, + HeatFluxBC, + TemperatureBC, + NaturalConvectionVerticalSpec, +) from tidy3d.components.tcad.generation_recombination import ( AugerRecombination, RadiativeRecombination, From 4e3013135e3ab74ff8fcdaafe872cfe9bd1d4c6c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Damiano=20Franz=C3=B2?= Date: Mon, 28 Jul 2025 19:18:43 +0200 Subject: [PATCH 3/9] Further improvements --- tidy3d/components/tcad/boundary/heat.py | 90 +++++++++++++++++++------ tidy3d/constants.py | 20 ++++++ 2 files changed, 89 insertions(+), 21 deletions(-) diff --git a/tidy3d/components/tcad/boundary/heat.py b/tidy3d/components/tcad/boundary/heat.py index 22e3567691..921f898031 100644 --- a/tidy3d/components/tcad/boundary/heat.py +++ b/tidy3d/components/tcad/boundary/heat.py @@ -2,11 +2,23 @@ from __future__ import annotations +from typing import Union + import pydantic.v1 as pd -from typing import Union from tidy3d.components.tcad.boundary.abstract import HeatChargeBC -from tidy3d.constants import HEAT_FLUX, HEAT_TRANSFER_COEFF, KELVIN +from tidy3d.constants import ( + ACCELERATION, + DENSITY, + DYNAMIC_VISCOSITY, + HEAT_FLUX, + HEAT_TRANSFER_COEFF, + KELVIN, + MICROMETER, + SPECIFIC_HEAT, + THERMAL_CONDUCTIVITY, + THERMAL_EXPANSIVITY, +) class TemperatureBC(HeatChargeBC): @@ -41,16 +53,6 @@ class HeatFluxBC(HeatChargeBC): ) -# TODO Should I adapt to Tidy3D units (e.g. micrometers), or keep the original ones? -THERMAL_CONDUCTIVITY_UNITS = "W/(m*K)" -DYNAMIC_VISCOSITY_UNITS = "Pa*s" -SPECIFIC_HEAT_UNITS = "J/(kg*K)" -DENSITY_UNITS = "kg/m**3" -THERMAL_EXPANSIVITY_UNITS = "1/K" -LENGTH_UNITS = "m" -ACCELERATION_UNITS = "m/s**2" - - class NaturalConvectionVerticalSpec(HeatChargeBC): """ Specification for natural convection from a vertical plate. @@ -65,41 +67,87 @@ class NaturalConvectionVerticalSpec(HeatChargeBC): fluid_k: pd.NonNegativeFloat = pd.Field( title="Fluid Thermal Conductivity", description="Thermal conductivity (k) of the fluid.", - units=THERMAL_CONDUCTIVITY_UNITS, + units=THERMAL_CONDUCTIVITY, ) fluid_mu: pd.NonNegativeFloat = pd.Field( title="Fluid Dynamic Viscosity", description="Dynamic viscosity (μ) of the fluid.", - units=DYNAMIC_VISCOSITY_UNITS, + units=DYNAMIC_VISCOSITY, ) fluid_Cp: pd.NonNegativeFloat = pd.Field( title="Fluid Specific Heat", - description="Specific heat capacity (Cp) of the fluid at constant pressure.", - units=SPECIFIC_HEAT_UNITS, + description="Specific heat of the fluid at constant pressure.", + units=SPECIFIC_HEAT, ) fluid_rho: pd.NonNegativeFloat = pd.Field( title="Fluid Density", description="Density (ρ) of the fluid.", - units=DENSITY_UNITS, + units=DENSITY, ) fluid_beta: pd.NonNegativeFloat = pd.Field( title="Fluid Thermal Expansivity", description="Thermal expansion coefficient (β) of the fluid.", - units=THERMAL_EXPANSIVITY_UNITS, + units=THERMAL_EXPANSIVITY, ) plate_L: pd.NonNegativeFloat = pd.Field( title="Plate Characteristic Length", description="Characteristic length (L), defined as the height of the vertical plate.", - units=LENGTH_UNITS, + units=MICROMETER, ) gravity: pd.NonNegativeFloat = pd.Field( - default=9.80665, + default=9.80665 * 1e6, title="Gravitational Acceleration", description="Gravitational acceleration (g).", - units=ACCELERATION_UNITS, + units=ACCELERATION, ) + def from_si_units( + fluid_k: pd.NonNegativeFloat, + fluid_mu: pd.NonNegativeFloat, + fluid_Cp: pd.NonNegativeFloat, + fluid_rho: pd.NonNegativeFloat, + fluid_beta: pd.NonNegativeFloat, + plate_L: pd.NonNegativeFloat, + gravity: float = 9.80665, + ): + """ + Create an instance from standard SI units. + + Args: + fluid_k: Thermal conductivity in [W/(m*K)]. + fluid_mu: Dynamic viscosity in [Pa*s]. + fluid_Cp: Specific heat in [J/(kg*K)]. + fluid_rho: Density in [kg/m**3]. + fluid_beta: Thermal expansivity in [1/K]. + plate_L: Plate characteristic length in [m]. + gravity: Gravitational acceleration in [m/s**2]. + + Returns: + An instance of NaturalConvectionVerticalSpec with all values + converted to Tidy3D's internal unit system. + """ + + # --- Apply conversion factors --- + # value_tidy = value_si * factor + k_tidy = fluid_k / 1e6 # W/(m*K) -> W/(um*K) + mu_tidy = fluid_mu / 1e6 # Pa*s -> kg/(um*s) + cp_tidy = fluid_Cp * 1e12 # J/(kg*K) -> um**2/(s**2*K) + rho_tidy = fluid_rho / 1e18 # kg/m**3 -> kg/um**3 + beta_tidy = fluid_beta # 1/K -> 1/K (no change) + L_tidy = plate_L * 1e6 # m -> um + g_tidy = gravity * 1e6 # m/s**2 -> um/s**2 + + return NaturalConvectionVerticalSpec( + fluid_k=k_tidy, + fluid_mu=mu_tidy, + fluid_Cp=cp_tidy, + fluid_rho=rho_tidy, + fluid_beta=beta_tidy, + plate_L=L_tidy, + gravity=g_tidy, + ) + def _compute_parameters(self): # Calculate the Rayleigh Number (Ra_L) rayleigh_numerator_notemp = ( diff --git a/tidy3d/constants.py b/tidy3d/constants.py index 7edcf22278..5d78c372d5 100644 --- a/tidy3d/constants.py +++ b/tidy3d/constants.py @@ -221,6 +221,26 @@ Amperes per square micrometer """ +DYNAMIC_VISCOSITY = "kg/(um*s)" +""" +Amperes per square micrometer +""" + +SPECIFIC_HEAT = "um^2/(s^2*K)" +""" +Square micrometers per (square second Kelvin). +""" + +THERMAL_EXPANSIVITY = "1/K" +""" +Inverse Kelvin. +""" + +ACCELERATION = "um/s^2" +""" +Acceleration unit. +""" + LARGE_NUMBER = 1e10 """ Large number used for comparing infinity. From fa13026003341764dea0f163a0f1ad95d946fdc5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Damiano=20Franz=C3=B2?= Date: Tue, 29 Jul 2025 16:42:35 +0200 Subject: [PATCH 4/9] Removing unnecessary update_forward_refs() --- tidy3d/components/tcad/boundary/heat.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/tidy3d/components/tcad/boundary/heat.py b/tidy3d/components/tcad/boundary/heat.py index 921f898031..04f58fd24d 100644 --- a/tidy3d/components/tcad/boundary/heat.py +++ b/tidy3d/components/tcad/boundary/heat.py @@ -204,8 +204,3 @@ class ConvectionBC(HeatChargeBC): description=f"Heat flux value in units of {HEAT_TRANSFER_COEFF}.", units=HEAT_TRANSFER_COEFF, ) - - -# TODO Maybe I should find a better place for these lines. -NaturalConvectionVerticalSpec.update_forward_refs() -ConvectionBC.update_forward_refs() From 65416202e81b66f5658375e4989aa437968b52a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Damiano=20Franz=C3=B2?= Date: Tue, 29 Jul 2025 17:43:38 +0200 Subject: [PATCH 5/9] Fix style with ruff format --- tidy3d/__init__.py | 4 ++-- tidy3d/components/tcad/types.py | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index f7bf678566..70b98986c5 100644 --- a/tidy3d/__init__.py +++ b/tidy3d/__init__.py @@ -24,6 +24,7 @@ from tidy3d.components.spice.sources.dc import DCCurrentSource, DCVoltageSource from tidy3d.components.spice.sources.types import VoltageSourceType from tidy3d.components.tcad.analysis.heat_simulation_type import UnsteadyHeatAnalysis, UnsteadySpec +from tidy3d.components.tcad.boundary.heat import NaturalConvectionVerticalSpec from tidy3d.components.tcad.boundary.specification import ( HeatBoundarySpec, HeatChargeBoundarySpec, @@ -68,7 +69,6 @@ CaugheyThomasMobility, ConstantMobilityModel, ConvectionBC, - NaturalConvectionVerticalSpec, CurrentBC, HeatFluxBC, HeatFromElectricSource, @@ -606,6 +606,7 @@ def set_logging_level(level: str) -> None: "ModulationSpec", "Monitor", "MultiPhysicsMedium", + "NaturalConvectionVerticalSpec", "NedeljkovicSorefMashanovich", "NonlinearModel", "NonlinearSpec", @@ -715,5 +716,4 @@ def set_logging_level(level: str) -> None: "set_logging_console", "set_logging_file", "wavelengths", - "NaturalConvectionVerticalSpec", ] diff --git a/tidy3d/components/tcad/types.py b/tidy3d/components/tcad/types.py index 38cb64c74b..5926eba088 100644 --- a/tidy3d/components/tcad/types.py +++ b/tidy3d/components/tcad/types.py @@ -8,7 +8,6 @@ ConvectionBC, HeatFluxBC, TemperatureBC, - NaturalConvectionVerticalSpec, ) from tidy3d.components.tcad.generation_recombination import ( AugerRecombination, From e7036eb427164970fd8761725984734ca9f1925b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Damiano=20Franz=C3=B2?= Date: Tue, 29 Jul 2025 17:46:51 +0200 Subject: [PATCH 6/9] Minor typo --- tidy3d/constants.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tidy3d/constants.py b/tidy3d/constants.py index 5d78c372d5..f637e451eb 100644 --- a/tidy3d/constants.py +++ b/tidy3d/constants.py @@ -223,7 +223,7 @@ DYNAMIC_VISCOSITY = "kg/(um*s)" """ -Amperes per square micrometer +Kilograms per (micrometer second) """ SPECIFIC_HEAT = "um^2/(s^2*K)" From c3ea537173dcfd973e54874ec8bf54d94b686266 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Damiano=20Franz=C3=B2?= Date: Sat, 2 Aug 2025 12:31:45 +0200 Subject: [PATCH 7/9] Update naming and minor fixes --- tidy3d/__init__.py | 5 ++--- tidy3d/components/tcad/boundary/heat.py | 13 +++++++------ 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index 70b98986c5..a7ec9f24db 100644 --- a/tidy3d/__init__.py +++ b/tidy3d/__init__.py @@ -24,7 +24,7 @@ from tidy3d.components.spice.sources.dc import DCCurrentSource, DCVoltageSource from tidy3d.components.spice.sources.types import VoltageSourceType from tidy3d.components.tcad.analysis.heat_simulation_type import UnsteadyHeatAnalysis, UnsteadySpec -from tidy3d.components.tcad.boundary.heat import NaturalConvectionVerticalSpec +from tidy3d.components.tcad.boundary.heat import VerticalNaturalConvectionCoeffModel from tidy3d.components.tcad.boundary.specification import ( HeatBoundarySpec, HeatChargeBoundarySpec, @@ -406,7 +406,6 @@ def set_logging_level(level: str) -> None: Transformed.update_forward_refs() ClipOperation.update_forward_refs() GeometryGroup.update_forward_refs() -ConvectionBC.update_forward_refs() __all__ = [ "C_0", @@ -606,7 +605,6 @@ def set_logging_level(level: str) -> None: "ModulationSpec", "Monitor", "MultiPhysicsMedium", - "NaturalConvectionVerticalSpec", "NedeljkovicSorefMashanovich", "NonlinearModel", "NonlinearSpec", @@ -696,6 +694,7 @@ def set_logging_level(level: str) -> None: "UnsteadyHeatAnalysis", "UnsteadySpec", "Updater", + "VerticalNaturalConvectionCoeffModel", "VisualizationSpec", "VoltageBC", "VoltageSourceType", diff --git a/tidy3d/components/tcad/boundary/heat.py b/tidy3d/components/tcad/boundary/heat.py index 04f58fd24d..84276bb27b 100644 --- a/tidy3d/components/tcad/boundary/heat.py +++ b/tidy3d/components/tcad/boundary/heat.py @@ -6,6 +6,7 @@ import pydantic.v1 as pd +from tidy3d.components.base import Tidy3dBaseModel from tidy3d.components.tcad.boundary.abstract import HeatChargeBC from tidy3d.constants import ( ACCELERATION, @@ -53,14 +54,14 @@ class HeatFluxBC(HeatChargeBC): ) -class NaturalConvectionVerticalSpec(HeatChargeBC): +class VerticalNaturalConvectionCoeffModel(Tidy3dBaseModel): """ Specification for natural convection from a vertical plate. This class calculates the heat transfer coefficient (h) based on fluid properties and an expected temperature difference, then provides these values as 'base' and 'exponent' for a generalized heat flux equation - q = base * (T_surf - T_fluid)^exponent. + q = base * (T_surf - T_fluid)^exponent + base * (T_surf - T_fluid). """ # --- Input Parameters --- @@ -124,7 +125,7 @@ def from_si_units( gravity: Gravitational acceleration in [m/s**2]. Returns: - An instance of NaturalConvectionVerticalSpec with all values + An instance of VerticalNaturalConvectionCoeffModel with all values converted to Tidy3D's internal unit system. """ @@ -138,7 +139,7 @@ def from_si_units( L_tidy = plate_L * 1e6 # m -> um g_tidy = gravity * 1e6 # m/s**2 -> um/s**2 - return NaturalConvectionVerticalSpec( + return VerticalNaturalConvectionCoeffModel( fluid_k=k_tidy, fluid_mu=mu_tidy, fluid_Cp=cp_tidy, @@ -175,7 +176,7 @@ def _compute_parameters(self): h_factor_non_linear = (0.387 * Ra_L ** (1 / 6)) / pr_denominator else: raise ValueError( - f"Ra_L={Ra_L} should be smaller than 1e13 for NaturalConvectionVerticalSpec" + f"Ra_L={Ra_L} should be smaller than 1e13 for VerticalNaturalConvectionCoeffModel" ) h = (self.fluid_k / self.plate_L) * h_factor_linear @@ -199,7 +200,7 @@ class ConvectionBC(HeatChargeBC): units=KELVIN, ) - transfer_coeff: Union[pd.NonNegativeFloat, NaturalConvectionVerticalSpec] = pd.Field( + transfer_coeff: Union[pd.NonNegativeFloat, VerticalNaturalConvectionCoeffModel] = pd.Field( title="Heat Transfer Coefficient", description=f"Heat flux value in units of {HEAT_TRANSFER_COEFF}.", units=HEAT_TRANSFER_COEFF, From 144650fb37c507eb711d47948c2959515778ee08 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Damiano=20Franz=C3=B2?= Date: Tue, 5 Aug 2025 17:06:22 +0200 Subject: [PATCH 8/9] Addressing last comments --- tidy3d/components/material/tcad/heat.py | 98 +++++++++++++++- tidy3d/components/tcad/boundary/heat.py | 108 +++--------------- .../components/tcad/simulation/heat_charge.py | 89 +++++++++++++++ tidy3d/constants.py | 5 + 4 files changed, 205 insertions(+), 95 deletions(-) diff --git a/tidy3d/components/material/tcad/heat.py b/tidy3d/components/material/tcad/heat.py index a33925f210..42f5fbf3a6 100644 --- a/tidy3d/components/material/tcad/heat.py +++ b/tidy3d/components/material/tcad/heat.py @@ -10,8 +10,11 @@ from tidy3d.components.base import Tidy3dBaseModel from tidy3d.constants import ( DENSITY, + DYNAMIC_VISCOSITY, + SPECIFIC_HEAT, SPECIFIC_HEAT_CAPACITY, THERMAL_CONDUCTIVITY, + THERMAL_EXPANSIVITY, ) @@ -46,11 +49,100 @@ class FluidMedium(AbstractHeatMedium): """Fluid medium. Heat simulations will not solve for temperature in a structure that has a medium with this 'heat_spec'. - Example - ------- - >>> solid = FluidMedium() + The full set of parameters is primarily intended for calculations involving natural + convection, where they are used to determine the heat transfer coefficient. + In the current version, these specific properties may not be utilized for + other boundary condition types. + + Attributes + ---------- + thermal_conductivity : float, optional + Thermal conductivity ($k$) of the fluid in $W/(\\mu m \\cdot K)$. + viscosity : float, optional + Dynamic viscosity ($\\mu$) of the fluid in $kg/(\\mu m \\cdot s)$. + specific_heat : float, optional + Specific heat ($c_p$) of the fluid in $\\mu m^2/(s^2 \\cdot K)$. + density : float, optional + Density ($\rho$) of the fluid in $kg/\\mu m^3$. + expansivity : float, optional + Thermal expansion coefficient ($\beta$) of the fluid in $1/K$. + + Examples + -------- + >>> # It is most convenient to define the fluid from standard SI units + >>> # using the `from_si_units` classmethod. + >>> # The following defines air at approximately 20°C. + >>> air_from_si = FluidMedium.from_si_units( + ... thermal_conductivity=0.0257, # Unit: W/(m*K) + ... viscosity=1.81e-5, # Unit: Pa*s + ... specific_heat=1005, # Unit: J/(kg*K) + ... density=1.204, # Unit: kg/m^3 + ... expansivity=1/293.15 # Unit: 1/K + ... ) + + >>> # One can also define the medium directly in Tidy3D units. + >>> # The following is equivalent to the example above. + >>> air_direct = FluidMedium( + ... thermal_conductivity=2.57e-8, + ... viscosity=1.81e-11, + ... specific_heat=1.005e+15, + ... density=1.204e-18, + ... expansivity=0.00341 + ... ) """ + thermal_conductivity: pd.NonNegativeFloat = pd.Field( + default=None, + title="Fluid Thermal Conductivity", + description="Thermal conductivity (k) of the fluid.", + units=THERMAL_CONDUCTIVITY, + ) + viscosity: pd.NonNegativeFloat = pd.Field( + default=None, + title="Fluid Dynamic Viscosity", + description="Dynamic viscosity (μ) of the fluid.", + units=DYNAMIC_VISCOSITY, + ) + specific_heat: pd.NonNegativeFloat = pd.Field( + default=None, + title="Fluid Specific Heat", + description="Specific heat of the fluid at constant pressure.", + units=SPECIFIC_HEAT, + ) + density: pd.NonNegativeFloat = pd.Field( + default=None, + title="Fluid Density", + description="Density (ρ) of the fluid.", + units=DENSITY, + ) + expansivity: pd.NonNegativeFloat = pd.Field( + default=None, + title="Fluid Thermal Expansivity", + description="Thermal expansion coefficient (β) of the fluid.", + units=THERMAL_EXPANSIVITY, + ) + + def from_si_units( + thermal_conductivity: pd.NonNegativeFloat, + viscosity: pd.NonNegativeFloat, + specific_heat: pd.NonNegativeFloat, + density: pd.NonNegativeFloat, + expansivity: pd.NonNegativeFloat, + ): + thermal_conductivity_tidy = thermal_conductivity / 1e6 # W/(m*K) -> W/(um*K) + viscosity_tidy = viscosity / 1e6 # Pa*s -> kg/(um*s) + specific_heat_tidy = specific_heat * 1e12 # J/(kg*K) -> um**2/(s**2*K) + density_tidy = density / 1e18 # kg/m**3 -> kg/um**3 + expansivity_tidy = expansivity # 1/K -> 1/K (no change) + + return FluidMedium( + thermal_conductivity=thermal_conductivity_tidy, + viscosity=viscosity_tidy, + specific_heat=specific_heat_tidy, + density=density_tidy, + expansivity=expansivity_tidy, + ) + class FluidSpec(FluidMedium): """Fluid medium class for backwards compatibility""" diff --git a/tidy3d/components/tcad/boundary/heat.py b/tidy3d/components/tcad/boundary/heat.py index 84276bb27b..3bc26c5165 100644 --- a/tidy3d/components/tcad/boundary/heat.py +++ b/tidy3d/components/tcad/boundary/heat.py @@ -7,18 +7,15 @@ import pydantic.v1 as pd from tidy3d.components.base import Tidy3dBaseModel +from tidy3d.components.material.tcad.heat import FluidMedium from tidy3d.components.tcad.boundary.abstract import HeatChargeBC from tidy3d.constants import ( ACCELERATION, - DENSITY, - DYNAMIC_VISCOSITY, + GRAV_ACC, HEAT_FLUX, HEAT_TRANSFER_COEFF, KELVIN, MICROMETER, - SPECIFIC_HEAT, - THERMAL_CONDUCTIVITY, - THERMAL_EXPANSIVITY, ) @@ -64,64 +61,35 @@ class VerticalNaturalConvectionCoeffModel(Tidy3dBaseModel): q = base * (T_surf - T_fluid)^exponent + base * (T_surf - T_fluid). """ - # --- Input Parameters --- - fluid_k: pd.NonNegativeFloat = pd.Field( - title="Fluid Thermal Conductivity", - description="Thermal conductivity (k) of the fluid.", - units=THERMAL_CONDUCTIVITY, + medium: FluidMedium = pd.Field( + default=None, + title="Interface medium", + description="Medium to use for heat transfer coefficient.", ) - fluid_mu: pd.NonNegativeFloat = pd.Field( - title="Fluid Dynamic Viscosity", - description="Dynamic viscosity (μ) of the fluid.", - units=DYNAMIC_VISCOSITY, - ) - fluid_Cp: pd.NonNegativeFloat = pd.Field( - title="Fluid Specific Heat", - description="Specific heat of the fluid at constant pressure.", - units=SPECIFIC_HEAT, - ) - fluid_rho: pd.NonNegativeFloat = pd.Field( - title="Fluid Density", - description="Density (ρ) of the fluid.", - units=DENSITY, - ) - fluid_beta: pd.NonNegativeFloat = pd.Field( - title="Fluid Thermal Expansivity", - description="Thermal expansion coefficient (β) of the fluid.", - units=THERMAL_EXPANSIVITY, - ) - plate_L: pd.NonNegativeFloat = pd.Field( + + plate_length: pd.NonNegativeFloat = pd.Field( title="Plate Characteristic Length", description="Characteristic length (L), defined as the height of the vertical plate.", units=MICROMETER, ) gravity: pd.NonNegativeFloat = pd.Field( - default=9.80665 * 1e6, + default=GRAV_ACC, title="Gravitational Acceleration", description="Gravitational acceleration (g).", units=ACCELERATION, ) def from_si_units( - fluid_k: pd.NonNegativeFloat, - fluid_mu: pd.NonNegativeFloat, - fluid_Cp: pd.NonNegativeFloat, - fluid_rho: pd.NonNegativeFloat, - fluid_beta: pd.NonNegativeFloat, - plate_L: pd.NonNegativeFloat, - gravity: float = 9.80665, + medium: FluidMedium, + plate_length: pd.NonNegativeFloat, + gravity: pd.NonNegativeFloat = GRAV_ACC * 1e-6, ): """ Create an instance from standard SI units. Args: - fluid_k: Thermal conductivity in [W/(m*K)]. - fluid_mu: Dynamic viscosity in [Pa*s]. - fluid_Cp: Specific heat in [J/(kg*K)]. - fluid_rho: Density in [kg/m**3]. - fluid_beta: Thermal expansivity in [1/K]. - plate_L: Plate characteristic length in [m]. + plate_length: Plate characteristic length in [m]. gravity: Gravitational acceleration in [m/s**2]. Returns: @@ -131,59 +99,15 @@ def from_si_units( # --- Apply conversion factors --- # value_tidy = value_si * factor - k_tidy = fluid_k / 1e6 # W/(m*K) -> W/(um*K) - mu_tidy = fluid_mu / 1e6 # Pa*s -> kg/(um*s) - cp_tidy = fluid_Cp * 1e12 # J/(kg*K) -> um**2/(s**2*K) - rho_tidy = fluid_rho / 1e18 # kg/m**3 -> kg/um**3 - beta_tidy = fluid_beta # 1/K -> 1/K (no change) - L_tidy = plate_L * 1e6 # m -> um + plate_length_tidy = plate_length * 1e6 # m -> um g_tidy = gravity * 1e6 # m/s**2 -> um/s**2 return VerticalNaturalConvectionCoeffModel( - fluid_k=k_tidy, - fluid_mu=mu_tidy, - fluid_Cp=cp_tidy, - fluid_rho=rho_tidy, - fluid_beta=beta_tidy, - plate_L=L_tidy, + medium=medium, + plate_length=plate_length_tidy, gravity=g_tidy, ) - def _compute_parameters(self): - # Calculate the Rayleigh Number (Ra_L) - rayleigh_numerator_notemp = ( - self.gravity * self.fluid_beta * self.fluid_rho**2 * self.fluid_Cp * self.plate_L**3 - ) - rayleigh_denominator = self.fluid_mu * self.fluid_k - Ra_L = rayleigh_numerator_notemp / rayleigh_denominator - - # Calculate the denominator term from the Nusselt number correlation - # This term is related to the Prandtl Number (Pr = mu * Cp / k) - pr_term_inner_num = 0.492 * self.fluid_k - pr_term_inner_den = self.fluid_mu * self.fluid_Cp - - pr_term_inner = pr_term_inner_num / pr_term_inner_den - pr_denominator = (1 + pr_term_inner ** (9 / 16)) ** (4 / 9) - - # Select formula based on flow regime (determined by Ra_L) - if Ra_L <= 1e9: - # Laminar Flow - h_factor_linear = 0.68 - h_factor_non_linear = (0.670 * Ra_L ** (1 / 6)) / pr_denominator - elif 1e9 < Ra_L <= 1e13: - # Turbulent Flow - h_factor_linear = 0.825 - h_factor_non_linear = (0.387 * Ra_L ** (1 / 6)) / pr_denominator - else: - raise ValueError( - f"Ra_L={Ra_L} should be smaller than 1e13 for VerticalNaturalConvectionCoeffModel" - ) - - h = (self.fluid_k / self.plate_L) * h_factor_linear - h_nonlinear = (self.fluid_k / self.plate_L) * h_factor_non_linear - exponent = 1 + 1 / 6 - return h, h_nonlinear, exponent - class ConvectionBC(HeatChargeBC): """Convective thermal boundary conditions. diff --git a/tidy3d/components/tcad/simulation/heat_charge.py b/tidy3d/components/tcad/simulation/heat_charge.py index 7b4b693c74..3631a8f241 100644 --- a/tidy3d/components/tcad/simulation/heat_charge.py +++ b/tidy3d/components/tcad/simulation/heat_charge.py @@ -9,6 +9,8 @@ import numpy as np import pydantic.v1 as pd +from tidy3d import FluidMedium, VerticalNaturalConvectionCoeffModel + try: from matplotlib import colormaps except ImportError: @@ -455,6 +457,93 @@ def check_voltage_array_if_capacitance(cls, values): ) return values + @pd.root_validator(skip_on_failure=True) + def check_natural_convection_bc(cls, values): + """Make sure that natural convection BCs are defined correctly.""" + boundary_spec = values.get("boundary_spec") + if not boundary_spec: + return values + + structures = values["structures"] + boundary_spec = values["boundary_spec"] + bg_medium = values["medium"] + + # Create mappings for easy lookup of media and structures by name. + media = {s.medium.name: s.medium for s in structures if s.medium.name} + if bg_medium and bg_medium.name: + media[bg_medium.name] = bg_medium + structures_map = {s.name: s for s in structures if s.name} + + def check_fluid_medium_attr(fluid_medium): + if ( + (fluid_medium.thermal_conductivity is None) + or (fluid_medium.viscosity is None) + or (fluid_medium.specific_heat is None) + or (fluid_medium.density is None) + or (fluid_medium.expansivity is None) + ): + raise SetupError( + f"Boundary spec at index {i}: The fluid medium at the natural convection interface " + f"must have 'thermal_conductivity', 'viscosity', 'specific_heat', 'density' and 'expansivity' defined." + ) + + for i, bc in enumerate(boundary_spec): + if not ( + isinstance(bc.condition, ConvectionBC) + and isinstance(bc.condition.transfer_coeff, VerticalNaturalConvectionCoeffModel) + ): + continue + + natural_conv_model = bc.condition.transfer_coeff + placement = bc.placement + + # Case 1: The fluid medium is inferred from the placement interface. + # We use direct dictionary access, assuming 'names_exist_bcs' validator has already run. + if natural_conv_model.medium is None: + if isinstance(placement, MediumMediumInterface): + if len(placement.mediums) != 2: + continue + med1 = media[placement.mediums[0]] + med2 = media[placement.mediums[1]] + elif isinstance(placement, StructureStructureInterface): + if len(placement.structures) != 2: + continue + med1 = structures_map[placement.structures[0]].medium + med2 = structures_map[placement.structures[1]].medium + else: + raise SetupError( + f"Boundary spec at index {i}: 'VerticalNaturalConvectionCoeffModel' with no medium specified requires " + f"the 'placement' to be of type 'MediumMediumInterface' or 'StructureStructureInterface', " + f"but got '{type(placement).__name__}'." + ) + specs = [ + med1.heat if isinstance(med1, MultiPhysicsMedium) else med1, + med2.heat if isinstance(med2, MultiPhysicsMedium) else med2, + ] + + # Check for a single fluid in the interface. + is_fluid = [isinstance(s, FluidMedium) for s in specs] + if is_fluid.count(True) != 1: + raise SetupError( + f"Boundary spec at index {i}: A natural convection boundary at an interface " + f"must be between exactly one solid and one fluid medium. " + f"Found types '{type(specs[0]).__name__}' and '{type(specs[1]).__name__}'." + ) + fluid_medium = specs[is_fluid.index(True)] + check_fluid_medium_attr(fluid_medium) + + # Case 2: The fluid medium IS specified directly in the convection model. + else: + fluid_medium = natural_conv_model.medium + if not isinstance(fluid_medium, FluidMedium): + raise SetupError( + f"Boundary spec at index {i}: The medium '{fluid_medium.name}' specified in " + f"'VerticalNaturalConvectionCoeffModel' must be a fluid, but it has a heat " + f"spec of type '{type(fluid_medium).__name__}'." + ) + check_fluid_medium_attr(fluid_medium) + return values + @pd.validator("size", always=True) def check_zero_dim_domain(cls, val, values): """Error if heat domain have zero dimensions.""" diff --git a/tidy3d/constants.py b/tidy3d/constants.py index f637e451eb..a9a0e3808d 100644 --- a/tidy3d/constants.py +++ b/tidy3d/constants.py @@ -53,6 +53,11 @@ Boltzmann constant [eV/K] """ +GRAV_ACC = 9.80665 * 1e6 +""" +Gravitational acceleration (g) [um/s^2].", +""" + # floating point precisions dp_eps = np.finfo(np.float64).eps """ From 3494f6d6aef930253893c3d88a123f69c17a1695 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Damiano=20Franz=C3=B2?= Date: Tue, 5 Aug 2025 20:50:21 +0200 Subject: [PATCH 9/9] Small fix --- tidy3d/components/tcad/boundary/heat.py | 2 +- tidy3d/components/tcad/simulation/heat_charge.py | 4 ---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/tidy3d/components/tcad/boundary/heat.py b/tidy3d/components/tcad/boundary/heat.py index 3bc26c5165..59337a9998 100644 --- a/tidy3d/components/tcad/boundary/heat.py +++ b/tidy3d/components/tcad/boundary/heat.py @@ -81,8 +81,8 @@ class VerticalNaturalConvectionCoeffModel(Tidy3dBaseModel): ) def from_si_units( - medium: FluidMedium, plate_length: pd.NonNegativeFloat, + medium: FluidMedium = None, gravity: pd.NonNegativeFloat = GRAV_ACC * 1e-6, ): """ diff --git a/tidy3d/components/tcad/simulation/heat_charge.py b/tidy3d/components/tcad/simulation/heat_charge.py index 3631a8f241..f68975de10 100644 --- a/tidy3d/components/tcad/simulation/heat_charge.py +++ b/tidy3d/components/tcad/simulation/heat_charge.py @@ -501,13 +501,9 @@ def check_fluid_medium_attr(fluid_medium): # We use direct dictionary access, assuming 'names_exist_bcs' validator has already run. if natural_conv_model.medium is None: if isinstance(placement, MediumMediumInterface): - if len(placement.mediums) != 2: - continue med1 = media[placement.mediums[0]] med2 = media[placement.mediums[1]] elif isinstance(placement, StructureStructureInterface): - if len(placement.structures) != 2: - continue med1 = structures_map[placement.structures[0]].medium med2 = structures_map[placement.structures[1]].medium else: