diff --git a/flow360/__init__.py b/flow360/__init__.py index e88aa8af1..aa897df90 100644 --- a/flow360/__init__.py +++ b/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", diff --git a/flow360/component/simulation/models/material.py b/flow360/component/simulation/models/material.py index c2aaa5faf..37561f423 100644 --- a/flow360/component/simulation/models/material.py +++ b/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,6 +208,10 @@ 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 ------- @@ -95,6 +219,20 @@ class Air(MaterialBase): ... 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: diff --git a/flow360/component/simulation/translator/solver_translator.py b/flow360/component/simulation/translator/solver_translator.py index 57a8040a9..97766eb9b 100644 --- a/flow360/component/simulation/translator/solver_translator.py +++ b/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