FEAT: Add natural convection BC

Author: damiano-flex

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",
 ]

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()

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,