John/generic gas

flow360/__init__.py

@@ -41,6 +41,11 @@
 )
 from flow360.component.simulation.models.material import (
     Air,
+    NASA9CoefficientSet,
+    NASA9Coefficients,
+    # Legacy aliases for backward compatibility
+    NASAPolynomialCoefficientSet,
+    NASAPolynomialCoefficients,
     SolidMaterial,
     Sutherland,
     Water,
@@ -266,6 +271,10 @@
     "Folder",
     "ForcePerArea",
     "Air",
+    "NASA9CoefficientSet",
+    "NASA9Coefficients",
+    "NASAPolynomialCoefficientSet",
+    "NASAPolynomialCoefficients",
     "Sutherland",
     "SolidMaterial",
     "Slice",

flow360/component/simulation/models/material.py

@@ -1,6 +1,6 @@
 """Material classes for the simulation framework."""
 
-from typing import Literal, Optional, Union
+from typing import List, Literal, Optional, Union
 
 import pydantic as pd
 from numpy import sqrt
@@ -28,6 +28,126 @@ class MaterialBase(Flow360BaseModel):
     name: str = pd.Field()
 
 
+class NASA9CoefficientSet(Flow360BaseModel):
+    """
+    Represents a set of 9 NASA polynomial coefficients for a specific temperature range.
+
+    The NASA 9-coefficient polynomial (McBride et al., 2002) computes thermodynamic
+    properties as:
+
+    cp/R = a0*T^-2 + a1*T^-1 + a2 + a3*T + a4*T^2 + a5*T^3 + a6*T^4
+
+    h/RT = -a0*T^-2 + a1*ln(T)/T + a2 + (a3/2)*T + (a4/3)*T^2 + (a5/4)*T^3 + (a6/5)*T^4 + a7/T
+
+    s/R = -(a0/2)*T^-2 - a1*T^-1 + a2*ln(T) + a3*T + (a4/2)*T^2 + (a5/3)*T^3 + (a6/4)*T^4 + a8
+
+    Coefficients:
+    - a0-a6: cp polynomial coefficients
+    - a7: enthalpy integration constant
+    - a8: entropy integration constant
+
+    Example
+    -------
+
+    >>> fl.NASA9CoefficientSet(
+    ...     temperature_range_min=200.0 * fl.u.K,
+    ...     temperature_range_max=1000.0 * fl.u.K,
+    ...     coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    ... )
+
+    ====
+    """
+
+    temperature_range_min: AbsoluteTemperatureType = pd.Field(
+        description="Minimum temperature for which this coefficient set is valid."
+    )
+    temperature_range_max: AbsoluteTemperatureType = pd.Field(
+        description="Maximum temperature for which this coefficient set is valid."
+    )
+    coefficients: List[float] = pd.Field(
+        description="Nine NASA polynomial coefficients [a0, a1, a2, a3, a4, a5, a6, a7, a8]. "
+        "a0-a6 are cp/R polynomial coefficients, a7 is the enthalpy integration constant, "
+        "and a8 is the entropy integration constant."
+    )
+
+    @pd.model_validator(mode="after")
+    def validate_coefficients(self):
+        """Validate that exactly 9 coefficients are provided."""
+        if len(self.coefficients) != 9:
+            raise ValueError(
+                f"NASA 9-coefficient polynomial requires exactly 9 coefficients, "
+                f"got {len(self.coefficients)}"
+            )
+        return self
+
+
+class NASA9Coefficients(Flow360BaseModel):
+    """
+    NASA 9-coefficient polynomial coefficients for computing temperature-dependent thermodynamic properties.
+
+    Supports 1-5 temperature ranges with continuous boundaries. Defaults to a single temperature range.
+
+    Example
+    -------
+
+    Single temperature range (default):
+
+    >>> fl.NASA9Coefficients(
+    ...     temperature_ranges=[
+    ...         fl.NASA9CoefficientSet(
+    ...             temperature_range_min=200.0 * fl.u.K,
+    ...             temperature_range_max=6000.0 * fl.u.K,
+    ...             coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    ...         )
+    ...     ]
+    ... )
+
+    Multiple temperature ranges:
+
+    >>> fl.NASA9Coefficients(
+    ...     temperature_ranges=[
+    ...         fl.NASA9CoefficientSet(
+    ...             temperature_range_min=200.0 * fl.u.K,
+    ...             temperature_range_max=1000.0 * fl.u.K,
+    ...             coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    ...         ),
+    ...         fl.NASA9CoefficientSet(
+    ...             temperature_range_min=1000.0 * fl.u.K,
+    ...             temperature_range_max=6000.0 * fl.u.K,
+    ...             coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    ...         )
+    ...     ]
+    ... )
+
+    ====
+    """
+
+    temperature_ranges: List[NASA9CoefficientSet] = pd.Field(
+        min_length=1,
+        max_length=5,
+        description="List of NASA 9-coefficient sets for different temperature ranges. "
+        "Must be ordered by increasing temperature and be continuous. Maximum 5 ranges supported."
+    )
+
+    @pd.model_validator(mode="after")
+    def validate_temperature_continuity(self):
+        """Validate that temperature ranges are continuous and non-overlapping."""
+        for i in range(len(self.temperature_ranges) - 1):
+            current_max = self.temperature_ranges[i].temperature_range_max
+            next_min = self.temperature_ranges[i + 1].temperature_range_min
+            if current_max != next_min:
+                raise ValueError(
+                    f"Temperature ranges must be continuous: range {i} max "
+                    f"({current_max}) must equal range {i+1} min ({next_min})"
+                )
+        return self
+
+
+# Legacy aliases for backward compatibility during transition
+NASAPolynomialCoefficientSet = NASA9CoefficientSet
+NASAPolynomialCoefficients = NASA9Coefficients
+
+
 class Sutherland(Flow360BaseModel):
     """
     Represents Sutherland's law for calculating dynamic viscosity.
@@ -88,13 +208,31 @@ class Air(MaterialBase):
     This sets specific material properties for air,
     including dynamic viscosity, specific heat ratio, gas constant, and Prandtl number.
 
+    The thermodynamic properties can be specified using NASA 9-coefficient polynomials
+    for temperature-dependent specific heats. By default, coefficients are set to
+    reproduce a constant gamma=1.4 (calorically perfect gas).
+
     Example
     -------
 
     >>> fl.Air(
     ...     dynamic_viscosity=1.063e-05 * fl.u.Pa * fl.u.s
     ... )
 
+    With custom NASA 9-coefficient polynomial:
+
+    >>> fl.Air(
+    ...     nasa_9_coefficients=fl.NASA9Coefficients(
+    ...         temperature_ranges=[
+    ...             fl.NASA9CoefficientSet(
+    ...                 temperature_range_min=200.0 * fl.u.K,
+    ...                 temperature_range_max=6000.0 * fl.u.K,
+    ...                 coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    ...             )
+    ...         ]
+    ...     )
+    ... )
+
     ====
     """
 
@@ -113,6 +251,26 @@ class Air(MaterialBase):
             "model with standard atmospheric conditions."
         ),
     )
+    nasa_9_coefficients: NASA9Coefficients = pd.Field(
+        default_factory=lambda: NASA9Coefficients(
+            temperature_ranges=[
+                NASA9CoefficientSet(
+                    temperature_range_min=200.0 * u.K,
+                    temperature_range_max=6000.0 * u.K,
+                    # For constant gamma=1.4: cp/R = gamma/(gamma-1) = 1.4/0.4 = 3.5
+                    # In NASA9 format, constant cp/R is the a2 coefficient (index 2)
+                    # All other coefficients (inverse T terms, positive T terms, integration constants) are zero
+                    coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+                ),
+            ]
+        ),
+        description=(
+            "NASA 9-coefficient polynomial coefficients for computing temperature-dependent "
+            "thermodynamic properties (cp, enthalpy, entropy). Defaults to a single temperature "
+            "range with coefficients that reproduce constant gamma=1.4 (calorically perfect gas). "
+            "For air with gamma=1.4: cp/R = 3.5 (stored in a2)."
+        ),
+    )
 
     @property
     def specific_heat_ratio(self) -> pd.PositiveFloat:

flow360/component/simulation/translator/solver_translator.py

@@ -12,7 +12,11 @@
     compute_udf_dimensionalization_factor,
 )
 from flow360.component.simulation.framework.entity_base import EntityList
-from flow360.component.simulation.models.material import Sutherland
+from flow360.component.simulation.models.material import (
+    Air,
+    NASA9Coefficients,
+    Sutherland,
+)
 from flow360.component.simulation.models.solver_numerics import NoneSolver
 from flow360.component.simulation.models.surface_models import (
     Freestream,
@@ -1483,6 +1487,90 @@ def calculate_monitor_semaphore_hash(params: SimulationParams):
     return hasher.hexdigest()
 
 
+def translate_nasa9_coefficients(
+    nasa_coeffs: NASA9Coefficients,
+    reference_temperature,
+):
+    """
+    Translate and non-dimensionalize NASA 9-coefficient polynomial coefficients.
+
+    The NASA 9-coefficient format (McBride et al., 2002):
+        cp/R = a0*T^-2 + a1*T^-1 + a2 + a3*T + a4*T^2 + a5*T^3 + a6*T^4
+        h/RT = -a0*T^-2 + a1*ln(T)/T + a2 + (a3/2)*T + (a4/3)*T^2 + (a5/4)*T^3 + (a6/5)*T^4 + a7/T
+        s/R = -(a0/2)*T^-2 - a1*T^-1 + a2*ln(T) + a3*T + (a4/2)*T^2 + (a5/3)*T^3 + (a6/4)*T^4 + a8
+
+    Coefficients: [a0, a1, a2, a3, a4, a5, a6, a7, a8]
+        - a0-a6: cp polynomial coefficients
+        - a7: enthalpy integration constant
+        - a8: entropy integration constant
+
+    In the C++ solver, the internal non-dimensional temperature is computed as:
+        T_internal = p / (R* * rho) where R* = 1/gamma_ref
+    This gives: T_internal = gamma_ref * T_dim / T_ref
+
+    Non-dimensionalization transformations (where t_scale = T_ref / gamma_ref):
+        - a0 (T^-2): a0_nd = a0 * (gamma_ref / T_ref)^2
+        - a1 (T^-1): a1_nd = a1 * (gamma_ref / T_ref)
+        - a2 (T^0):  a2_nd = a2 (unchanged)
+        - a3 (T^1):  a3_nd = a3 * t_scale
+        - a4 (T^2):  a4_nd = a4 * t_scale^2
+        - a5 (T^3):  a5_nd = a5 * t_scale^3
+        - a6 (T^4):  a6_nd = a6 * t_scale^4
+        - a7 (1/T):  a7_nd = a7 * (gamma_ref / T_ref)
+        - a8 (const): a8_nd = a8 (unchanged)
+
+    Parameters
+    ----------
+    nasa_coeffs : NASA9Coefficients
+        NASA 9-coefficient polynomial coefficients with temperature ranges
+    reference_temperature : float
+        Reference temperature for non-dimensionalization (in K)
+
+    Returns
+    -------
+    dict
+        Non-dimensionalized NASA 9-coefficient data for Flow360.json
+    """
+
+    # Reference gamma used in Flow360 non-dimensionalization
+    gamma_ref = 1.4
+
+    # Scaling factors
+    t_scale = reference_temperature / gamma_ref  # For positive powers of T
+    t_scale_inv = gamma_ref / reference_temperature  # For negative powers of T
+
+    temperature_ranges = []
+    for coeff_set in nasa_coeffs.temperature_ranges:
+        # Temperature ranges need to be scaled by gamma_ref (T_internal = gamma_ref * T_dim / T_ref)
+        t_min_nd = gamma_ref * coeff_set.temperature_range_min.to("K").v.item() / reference_temperature
+        t_max_nd = gamma_ref * coeff_set.temperature_range_max.to("K").v.item() / reference_temperature
+
+        coeffs = coeff_set.coefficients
+
+        # Transform coefficients for T_internal = gamma*T_dim/T_ref
+        coeffs_nd = [
+            coeffs[0] * t_scale_inv**2,  # a0 (T^-2): scale by (gamma/T_ref)^2
+            coeffs[1] * t_scale_inv,  # a1 (T^-1): scale by (gamma/T_ref)
+            coeffs[2],  # a2 (constant): no scaling
+            coeffs[3] * t_scale,  # a3 (T^1): scale by (T_ref/gamma)
+            coeffs[4] * t_scale**2,  # a4 (T^2): scale by (T_ref/gamma)^2
+            coeffs[5] * t_scale**3,  # a5 (T^3): scale by (T_ref/gamma)^3
+            coeffs[6] * t_scale**4,  # a6 (T^4): scale by (T_ref/gamma)^4
+            coeffs[7] * t_scale_inv,  # a7 (enthalpy, 1/T): scale by (gamma/T_ref)
+            coeffs[8],  # a8 (entropy, constant): no scaling
+        ]
+
+        temperature_ranges.append(
+            {
+                "temperatureRangeMin": t_min_nd,
+                "temperatureRangeMax": t_max_nd,
+                "coefficients": coeffs_nd,
+            }
+        )
+
+    return {"temperatureRanges": temperature_ranges}
+
+
 # pylint: disable=too-many-statements
 # pylint: disable=too-many-branches
 # pylint: disable=too-many-locals
@@ -1534,6 +1622,16 @@ def get_solver_json(
             else op.material.dynamic_viscosity.in_base(input_params.flow360_unit_system).v.item()
         ),
     }
+
+    ##:: Step 2.5: Get NASA 9-coefficient polynomial for gas model (if Air material)
+    if not isinstance(op, LiquidOperatingCondition) and isinstance(op.thermal_state.material, Air):
+        # Get reference temperature for non-dimensionalization (freestream temperature)
+        reference_temperature = op.thermal_state.temperature.to("K").v.item()
+
+        translated["thermallyPerfectGasModel"] = translate_nasa9_coefficients(
+            op.thermal_state.material.nasa_9_coefficients,
+            reference_temperature,
+        )
     if (
         "reference_velocity_magnitude" in op.__class__.model_fields.keys()
         and op.reference_velocity_magnitude