diff --git a/flow360/__init__.py b/flow360/__init__.py index e88aa8af1..ee85c2e09 100644 --- a/flow360/__init__.py +++ b/flow360/__init__.py @@ -41,8 +41,15 @@ ) from flow360.component.simulation.models.material import ( Air, + FrozenSpecies, + NASA9CoefficientSet, + NASA9Coefficients, + # Legacy aliases for backward compatibility + NASAPolynomialCoefficientSet, + NASAPolynomialCoefficients, SolidMaterial, Sutherland, + ThermallyPerfectGas, Water, ) from flow360.component.simulation.models.solver_numerics import ( @@ -266,7 +273,13 @@ "Folder", "ForcePerArea", "Air", + "NASA9CoefficientSet", + "NASA9Coefficients", + "NASAPolynomialCoefficientSet", + "NASAPolynomialCoefficients", + "FrozenSpecies", "Sutherland", + "ThermallyPerfectGas", "SolidMaterial", "Slice", "Isosurface", diff --git a/flow360/component/simulation/models/material.py b/flow360/component/simulation/models/material.py index c2aaa5faf..b5d12dd54 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,227 @@ 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 FrozenSpecies(Flow360BaseModel): + """ + Represents a single gas species with NASA 9-coefficient thermodynamic properties. + + Used within :class:`ThermallyPerfectGas` to define multi-species gas mixtures + where each species contributes to the mixture properties weighted by mass fraction. + The term "frozen" indicates fixed mass fractions (non-reacting flow). + + Example + ------- + + >>> fl.FrozenSpecies( + ... name="N2", + ... 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] + ... ) + ... ] + ... ), + ... mass_fraction=0.7555 + ... ) + + ==== + """ + + name: str = pd.Field(description="Species name (e.g., 'N2', 'O2', 'Ar')") + nasa_9_coefficients: NASA9Coefficients = pd.Field( + description="NASA 9-coefficient polynomial for this species" + ) + mass_fraction: pd.PositiveFloat = pd.Field( + description="Mass fraction of this species (must sum to 1 across all species in mixture)" + ) + + +class ThermallyPerfectGas(Flow360BaseModel): + """ + Multi-species thermally perfect gas model. + + Combines NASA 9-coefficient polynomials from multiple species weighted by mass fraction. + All species must use the same temperature range boundaries. The mixture properties + are computed as mass-fraction-weighted averages of individual species properties. + + This model supports temperature-dependent specific heats (cp) while maintaining + fixed mass fractions (non-reacting flow). + + Example + ------- + + >>> fl.ThermallyPerfectGas( + ... species=[ + ... fl.FrozenSpecies(name="N2", nasa_9_coefficients=..., mass_fraction=0.7555), + ... fl.FrozenSpecies(name="O2", nasa_9_coefficients=..., mass_fraction=0.2316), + ... fl.FrozenSpecies(name="Ar", nasa_9_coefficients=..., mass_fraction=0.0129), + ... ] + ... ) + + ==== + """ + + species: List[FrozenSpecies] = pd.Field( + min_length=1, + description="List of species with their NASA 9 coefficients and mass fractions. " + "Mass fractions must sum to 1.0." + ) + + @pd.model_validator(mode="after") + def validate_mass_fractions_sum_to_one(self): + """Validate that mass fractions sum to 1.""" + total = sum(s.mass_fraction for s in self.species) + if not (0.999 <= total <= 1.001): # Allow small tolerance for floating point + raise ValueError(f"Mass fractions must sum to 1.0, got {total}") + return self + + @pd.model_validator(mode="after") + def validate_temperature_ranges_match(self): + """Validate all species have matching temperature range boundaries.""" + if len(self.species) < 2: + return self + ref_ranges = self.species[0].nasa_9_coefficients.temperature_ranges + for species in self.species[1:]: + ranges = species.nasa_9_coefficients.temperature_ranges + if len(ranges) != len(ref_ranges): + raise ValueError( + f"Species '{species.name}' has {len(ranges)} temperature ranges, " + f"but '{self.species[0].name}' has {len(ref_ranges)}. " + "All species must have the same number of temperature ranges." + ) + for i, (r1, r2) in enumerate(zip(ref_ranges, ranges)): + if r1.temperature_range_min != r2.temperature_range_min or \ + r1.temperature_range_max != r2.temperature_range_max: + raise ValueError( + f"Temperature range {i} boundaries mismatch between species " + f"'{self.species[0].name}' and '{species.name}'. " + "All species must use the same temperature range boundaries." + ) + return self + + class Sutherland(Flow360BaseModel): """ Represents Sutherland's law for calculating dynamic viscosity. @@ -88,6 +309,13 @@ 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). + + For multi-species gas mixtures, use the `thermally_perfect_gas` parameter which + combines species properties weighted by mass fraction. + Example ------- @@ -95,6 +323,32 @@ 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] + ... ) + ... ] + ... ) + ... ) + + With multi-species thermally perfect gas: + + >>> fl.Air( + ... thermally_perfect_gas=fl.ThermallyPerfectGas( + ... species=[ + ... fl.FrozenSpecies(name="N2", nasa_9_coefficients=..., mass_fraction=0.7555), + ... fl.FrozenSpecies(name="O2", nasa_9_coefficients=..., mass_fraction=0.2316), + ... fl.FrozenSpecies(name="Ar", nasa_9_coefficients=..., mass_fraction=0.0129), + ... ] + ... ) + ... ) + ==== """ @@ -113,18 +367,117 @@ 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). " + "Note: If thermally_perfect_gas is specified, it takes precedence over this field." + ), + ) + thermally_perfect_gas: Optional[ThermallyPerfectGas] = pd.Field( + default=None, + description=( + "Multi-species thermally perfect gas model. When specified, this takes precedence " + "over nasa_9_coefficients. Use this to define gas mixtures with multiple species, " + "each with their own NASA 9-coefficient polynomials and mass fractions. " + "The mixture properties are computed as mass-fraction-weighted averages." + ), + ) + + def get_specific_heat_ratio(self, temperature: AbsoluteTemperatureType) -> pd.PositiveFloat: + """ + Computes the specific heat ratio (gamma) at a given temperature from NASA polynomial. + + For thermally perfect gas, gamma = cp/cv = (cp/R) / (cp/R - 1) varies with temperature. + The cp/R is computed from the NASA 9-coefficient polynomial: + cp/R = a0*T^-2 + a1*T^-1 + a2 + a3*T + a4*T^2 + a5*T^3 + a6*T^4 + + Parameters + ---------- + temperature : AbsoluteTemperatureType + The temperature at which to compute gamma. + + Returns + ------- + pd.PositiveFloat + The specific heat ratio at the given temperature. + """ + T = temperature.to("K").v.item() + + # Get coefficients from the appropriate source + if self.thermally_perfect_gas is not None: + # For multi-species, combine coefficients by mass fraction + coeffs = [0.0] * 9 + for species in self.thermally_perfect_gas.species: + # Find the temperature range that contains T + for coeff_set in species.nasa_9_coefficients.temperature_ranges: + t_min = coeff_set.temperature_range_min.to("K").v.item() + t_max = coeff_set.temperature_range_max.to("K").v.item() + if t_min <= T <= t_max: + for i in range(9): + coeffs[i] += species.mass_fraction * coeff_set.coefficients[i] + break + else: + # Single-species: find the temperature range that contains T + coeffs = None + for coeff_set in self.nasa_9_coefficients.temperature_ranges: + t_min = coeff_set.temperature_range_min.to("K").v.item() + t_max = coeff_set.temperature_range_max.to("K").v.item() + if t_min <= T <= t_max: + coeffs = list(coeff_set.coefficients) + break + if coeffs is None: + # Fallback to first range if T is out of bounds + coeffs = list(self.nasa_9_coefficients.temperature_ranges[0].coefficients) + + # Compute cp/R from the polynomial + # cp/R = a0*T^-2 + a1*T^-1 + a2 + a3*T + a4*T^2 + a5*T^3 + a6*T^4 + cp_over_R = ( + coeffs[0] * T ** (-2) + + coeffs[1] * T ** (-1) + + coeffs[2] + + coeffs[3] * T + + coeffs[4] * T**2 + + coeffs[5] * T**3 + + coeffs[6] * T**4 + ) + + # cv/R = cp/R - 1 (for ideal gas: cp - cv = R) + cv_over_R = cp_over_R - 1 + + # gamma = cp/cv + gamma = cp_over_R / cv_over_R + + return gamma @property def specific_heat_ratio(self) -> pd.PositiveFloat: """ - Returns the specific heat ratio (gamma) for air. + Returns the specific heat ratio (gamma) for air at a reference temperature (298.15 K). + + For temperature-dependent gamma, use `get_specific_heat_ratio(temperature)` instead. Returns ------- pd.PositiveFloat - The specific heat ratio, typically 1.4 for air. + The specific heat ratio at 298.15 K. """ - return 1.4 + # Compute gamma at reference temperature (298.15 K) + return self.get_specific_heat_ratio(298.15 * u.K) @property def gas_constant(self) -> SpecificHeatCapacityType.Positive: @@ -179,6 +532,8 @@ def get_speed_of_sound(self, temperature: AbsoluteTemperatureType) -> VelocityTy """ Calculates the speed of sound in air at a given temperature. + For thermally perfect gas, uses the temperature-dependent gamma from the NASA polynomial. + Parameters ---------- temperature : AbsoluteTemperatureType @@ -190,7 +545,8 @@ def get_speed_of_sound(self, temperature: AbsoluteTemperatureType) -> VelocityTy The speed of sound at the specified temperature. """ temperature = temperature.to("K") - return sqrt(self.specific_heat_ratio * self.gas_constant * temperature) + gamma = self.get_specific_heat_ratio(temperature) + return sqrt(gamma * self.gas_constant * temperature) @pd.validate_call def get_dynamic_viscosity( diff --git a/flow360/component/simulation/simulation_params.py b/flow360/component/simulation/simulation_params.py index 5309d75af..2194b6ee8 100644 --- a/flow360/component/simulation/simulation_params.py +++ b/flow360/component/simulation/simulation_params.py @@ -107,6 +107,7 @@ _check_numerical_dissipation_factor_output, _check_parent_volume_is_rotating, _check_time_average_output, + _check_tpg_not_with_isentropic_solver, _check_unsteadiness_to_use_hybrid_model, _check_valid_models_for_liquid, ) @@ -544,6 +545,11 @@ def check_low_mach_preconditioner_output(self): """Only allow lowMachPreconditioner output field when the lowMachPreconditioner is enabled in the NS solver""" return _check_low_mach_preconditioner_output(self) + @pd.model_validator(mode="after") + def check_tpg_not_with_isentropic_solver(self): + """ThermallyPerfectGas model is not supported with the CompressibleIsentropic (4x4) solver""" + return _check_tpg_not_with_isentropic_solver(self) + @pd.model_validator(mode="after") @context_validator(context=CASE) def check_complete_boundary_condition_and_unknown_surface(self): diff --git a/flow360/component/simulation/translator/solver_translator.py b/flow360/component/simulation/translator/solver_translator.py index 57a8040a9..e217d95cd 100644 --- a/flow360/component/simulation/translator/solver_translator.py +++ b/flow360/component/simulation/translator/solver_translator.py @@ -12,7 +12,12 @@ 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, + ThermallyPerfectGas, +) from flow360.component.simulation.models.solver_numerics import NoneSolver from flow360.component.simulation.models.surface_models import ( Freestream, @@ -1483,6 +1488,292 @@ def calculate_monitor_semaphore_hash(params: SimulationParams): return hasher.hexdigest() +def _compute_a7_correction(coeffs: list, reference_temperature: float) -> float: + """ + Compute the corrected a7 coefficient for Flow360 non-dimensionalization. + + The NASA polynomial a7 (enthalpy integration constant) is calibrated for absolute + dimensional enthalpy (h=0 at some reference temperature like 0K or 298.15K). + However, Flow360's non-dimensionalization expects internal energy to be consistent + with the ideal gas equation of state: e = cv * T at T_ref. + + This function computes a corrected a7 so that at T_nd = 1 (T = T_ref): + h_nd/R = gamma/(gamma-1) (ideal gas non-dimensional enthalpy) + e_nd/R = 1/(gamma-1) (ideal gas non-dimensional internal energy) + + The gamma at T_ref is computed directly from the polynomial: + cp/R = a0*T^-2 + a1*T^-1 + a2 + a3*T + a4*T^2 + a5*T^3 + a6*T^4 + cv/R = cp/R - 1 + gamma = cp/cv = (cp/R) / (cp/R - 1) + + IMPORTANT: The ln(T) term requires special handling. Under non-dimensionalization + T -> T_nd = T/T_ref, we have ln(T) = ln(T_nd) + ln(T_ref). At T_nd = 1, ln(T_nd) = 0, + so the a1*ln(T_ref) contribution must be absorbed into a7 to maintain consistency: + a7_corrected = a7_for_h_target + a1*ln(T_ref) + + Parameters + ---------- + coeffs : list + Dimensional NASA 9 coefficients [a0, a1, a2, a3, a4, a5, a6, a7, a8] + reference_temperature : float + Reference temperature for non-dimensionalization (in K) + + Returns + ------- + float + Corrected dimensional a7 coefficient (with ln(T_ref) term absorbed) + """ + import math + + T = reference_temperature + + # Compute cp/R at T_ref from the polynomial + # cp/R = a0*T^-2 + a1*T^-1 + a2 + a3*T + a4*T^2 + a5*T^3 + a6*T^4 + cp_over_R = ( + coeffs[0] * T ** (-2) + + coeffs[1] * T ** (-1) + + coeffs[2] + + coeffs[3] * T + + coeffs[4] * T**2 + + coeffs[5] * T**3 + + coeffs[6] * T**4 + ) + + # cv/R = cp/R - 1 (for ideal gas: cp - cv = R) + cv_over_R = cp_over_R - 1 + + # gamma = cp/cv at T_ref + gamma = cp_over_R / cv_over_R + + # Compute h/R at T_ref without the a7 term (dimensional) + # h/R = -a0/T + a1*ln(T) + a2*T + (a3/2)*T^2 + (a4/3)*T^3 + (a5/4)*T^4 + (a6/5)*T^5 + h_no_a7 = ( + -coeffs[0] / T + + coeffs[1] * math.log(T) + + coeffs[2] * T + + (coeffs[3] / 2) * T**2 + + (coeffs[4] / 3) * T**3 + + (coeffs[5] / 4) * T**4 + + (coeffs[6] / 5) * T**5 + ) + + # Target h/R at T_ref for ideal gas with computed gamma + # h/R = cp/R * T = (gamma/(gamma-1)) * T + h_target = (gamma / (gamma - 1)) * T + + # Base a7 correction ensures h_polynomial = h_target at T_ref (dimensional) + a7_base = h_target - h_no_a7 + + # Absorb the a1*ln(T_ref) term into a7 + # Under non-dimensionalization, ln(T) = ln(T_nd * T_ref) = ln(T_nd) + ln(T_ref) + # At T_nd = 1, ln(T_nd) = 0, so the a1*ln(T_ref) constant must be in a7 + a7_corrected = a7_base + coeffs[1] * math.log(T) + + return a7_corrected + + +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/R = -a0/T + a1*ln(T) + a2*T + (a3/2)*T^2 + (a4/3)*T^3 + (a5/4)*T^4 + (a6/5)*T^5 + a7 + 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 + + IMPORTANT: The a7 coefficient is automatically corrected to ensure internal energy + consistency with Flow360's non-dimensionalization. The NASA polynomial a7 is + calibrated for absolute enthalpy, but Flow360 expects e = cv*T at T_ref. + + In the C++ solver, the internal non-dimensional temperature is computed as: + T_internal = p* * gamma / rho* where p* = p/(rho_ref*a_ref^2), rho* = rho/rho_ref + Using p = rho*R*T and a_ref^2 = gamma*R*T_ref, this gives: + T_internal = T_dim / T_ref (NOT gamma * T_dim / T_ref) + + Non-dimensionalization transformations (where t_scale = T_ref): + - a0 (T^-2): a0_nd = a0 / T_ref^2 + - a1 (T^-1): a1_nd = a1 / T_ref + - a2 (T^0): a2_nd = a2 (unchanged) + - a3 (T^1): a3_nd = a3 * T_ref + - a4 (T^2): a4_nd = a4 * T_ref^2 + - a5 (T^3): a5_nd = a5 * T_ref^3 + - a6 (T^4): a6_nd = a6 * T_ref^4 + - a7 (const): a7_nd = a7 / T_ref (after correction) + - 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 + """ + + # Scaling factors for non-dimensionalization + # T_internal = T_dim / T_ref, so t_scale = T_ref + t_scale = reference_temperature # For positive powers of T + t_scale_inv = 1.0 / reference_temperature # For negative powers of T + + temperature_ranges = [] + for coeff_set in nasa_coeffs.temperature_ranges: + # Temperature ranges: T_internal = T_dim / T_ref + t_min_nd = coeff_set.temperature_range_min.to("K").v.item() / reference_temperature + t_max_nd = coeff_set.temperature_range_max.to("K").v.item() / reference_temperature + + coeffs = list(coeff_set.coefficients) + + # Correct a7 for Flow360 non-dimensionalization + # This ensures internal energy consistency with e = cv*T at T_ref + coeffs[7] = _compute_a7_correction(coeffs, reference_temperature) + + # Transform coefficients for T_internal = T_dim / T_ref + coeffs_nd = [ + coeffs[0] * t_scale_inv**2, # a0 (T^-2): scale by (1/T_ref)^2 + coeffs[1] * t_scale_inv, # a1 (T^-1): scale by (1/T_ref) + coeffs[2], # a2 (constant): no scaling + coeffs[3] * t_scale, # a3 (T^1): scale by T_ref + coeffs[4] * t_scale**2, # a4 (T^2): scale by T_ref^2 + coeffs[5] * t_scale**3, # a5 (T^3): scale by T_ref^3 + coeffs[6] * t_scale**4, # a6 (T^4): scale by T_ref^4 + coeffs[7] * t_scale_inv, # a7 (enthalpy const): scale by (1/T_ref) + coeffs[8], # a8 (entropy, constant): no scaling + ] + + temperature_ranges.append( + { + "temperatureRangeMin": t_min_nd, + "temperatureRangeMax": t_max_nd, + "coefficients": coeffs_nd, + } + ) + + # Compute gasConstant = 1/gamma at T_nd=1 from the first temperature range + # This ensures correct non-dimensionalization for temperature output + first_coeffs = temperature_ranges[0]["coefficients"] + T_nd = 1.0 + cp_over_R = ( + first_coeffs[0] * T_nd ** (-2) + + first_coeffs[1] * T_nd ** (-1) + + first_coeffs[2] + + first_coeffs[3] * T_nd + + first_coeffs[4] * T_nd**2 + + first_coeffs[5] * T_nd**3 + + first_coeffs[6] * T_nd**4 + ) + gamma_at_Tref = cp_over_R / (cp_over_R - 1.0) + gas_constant = 1.0 / gamma_at_Tref + + return {"temperatureRanges": temperature_ranges, "gasConstant": gas_constant} + + +def translate_thermally_perfect_gas( + tpg: ThermallyPerfectGas, + reference_temperature: float, +): + """ + Translate multi-species thermally perfect gas to Flow360 format. + + Combines NASA 9-coefficient polynomials from multiple species weighted by mass fraction. + For each coefficient, the mixture value is computed as: + a_mixture[i] = sum(mass_fraction[j] * a[j][i] for all species j) + + This assumes all species share the same temperature range boundaries (validated in + ThermallyPerfectGas class). + + IMPORTANT: The a7 coefficient is automatically corrected to ensure internal energy + consistency with Flow360's non-dimensionalization. The NASA polynomial a7 is + calibrated for absolute enthalpy, but Flow360 expects e = cv*T at T_ref. + + Parameters + ---------- + tpg : ThermallyPerfectGas + Multi-species thermally perfect gas model with species and mass fractions + reference_temperature : float + Reference temperature for non-dimensionalization (in K) + + Returns + ------- + dict + Non-dimensionalized NASA 9-coefficient data for Flow360.json with mass-fraction + weighted coefficients + """ + # Use first species as template for temperature ranges (all validated to match) + ref_ranges = tpg.species[0].nasa_9_coefficients.temperature_ranges + + # Scaling factors for non-dimensionalization + t_scale = reference_temperature + t_scale_inv = 1.0 / reference_temperature + + temperature_ranges = [] + for range_idx, ref_range in enumerate(ref_ranges): + # Combine coefficients weighted by mass fraction + combined_coeffs = [0.0] * 9 + for species in tpg.species: + coeffs = species.nasa_9_coefficients.temperature_ranges[range_idx].coefficients + for i in range(9): + combined_coeffs[i] += species.mass_fraction * coeffs[i] + + # Correct a7 for Flow360 non-dimensionalization + # This ensures internal energy consistency with e = cv*T at T_ref + combined_coeffs[7] = _compute_a7_correction(combined_coeffs, reference_temperature) + + # Non-dimensionalize (same as translate_nasa9_coefficients) + coeffs_nd = [ + combined_coeffs[0] * t_scale_inv**2, # a0 (T^-2): scale by (1/T_ref)^2 + combined_coeffs[1] * t_scale_inv, # a1 (T^-1): scale by (1/T_ref) + combined_coeffs[2], # a2 (constant): no scaling + combined_coeffs[3] * t_scale, # a3 (T^1): scale by T_ref + combined_coeffs[4] * t_scale**2, # a4 (T^2): scale by T_ref^2 + combined_coeffs[5] * t_scale**3, # a5 (T^3): scale by T_ref^3 + combined_coeffs[6] * t_scale**4, # a6 (T^4): scale by T_ref^4 + combined_coeffs[7] * t_scale_inv, # a7 (enthalpy const): scale by (1/T_ref) + combined_coeffs[8], # a8 (entropy, constant): no scaling + ] + + # Temperature ranges: T_internal = T_dim / T_ref + t_min_nd = ref_range.temperature_range_min.to("K").v.item() / reference_temperature + t_max_nd = ref_range.temperature_range_max.to("K").v.item() / reference_temperature + + temperature_ranges.append( + { + "temperatureRangeMin": t_min_nd, + "temperatureRangeMax": t_max_nd, + "coefficients": coeffs_nd, + } + ) + + # Compute gasConstant = 1/gamma at T_nd=1 from the first temperature range + # This ensures correct non-dimensionalization for temperature output + first_coeffs = temperature_ranges[0]["coefficients"] + T_nd = 1.0 + cp_over_R = ( + first_coeffs[0] * T_nd ** (-2) + + first_coeffs[1] * T_nd ** (-1) + + first_coeffs[2] + + first_coeffs[3] * T_nd + + first_coeffs[4] * T_nd**2 + + first_coeffs[5] * T_nd**3 + + first_coeffs[6] * T_nd**4 + ) + gamma_at_Tref = cp_over_R / (cp_over_R - 1.0) + gas_constant = 1.0 / gamma_at_Tref + + return {"temperatureRanges": temperature_ranges, "gasConstant": gas_constant} + + # pylint: disable=too-many-statements # pylint: disable=too-many-branches # pylint: disable=too-many-locals @@ -1534,6 +1825,23 @@ 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() + + # Check for multi-species thermally perfect gas first, then fall back to single-species + if op.thermal_state.material.thermally_perfect_gas is not None: + translated["thermallyPerfectGasModel"] = translate_thermally_perfect_gas( + op.thermal_state.material.thermally_perfect_gas, + reference_temperature, + ) + else: + 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 diff --git a/flow360/component/simulation/validation/validation_simulation_params.py b/flow360/component/simulation/validation/validation_simulation_params.py index 76871be3b..d9ab48421 100644 --- a/flow360/component/simulation/validation/validation_simulation_params.py +++ b/flow360/component/simulation/validation/validation_simulation_params.py @@ -5,6 +5,7 @@ from typing import Type, Union, get_args from flow360.component.simulation.meshing_param.volume_params import CustomZones +from flow360.component.simulation.models.material import Air from flow360.component.simulation.models.solver_numerics import NoneSolver from flow360.component.simulation.models.surface_models import ( Inflow, @@ -598,3 +599,44 @@ def _check_actuator_disk_names(models): _check_actuator_disk_names(models) return models + + +def _check_tpg_not_with_isentropic_solver(params): + """ + Validate that ThermallyPerfectGas model is not used with CompressibleIsentropic solver. + + The CompressibleIsentropic solver (4x4 system) does not support thermally perfect gas + models for performance reasons. Users must use the full Compressible solver when using + temperature-dependent gas properties. + """ + # Check if CompressibleIsentropic solver is being used + uses_isentropic = False + if params.models: + for model in params.models: + if isinstance(model, Fluid): + if model.navier_stokes_solver.type_name == "CompressibleIsentropic": + uses_isentropic = True + break + + if not uses_isentropic: + return params + + # Check if TPG model is being used (via operating condition) + uses_tpg = False + if params.operating_condition is not None: + op = params.operating_condition + if hasattr(op, "thermal_state") and op.thermal_state is not None: + material = op.thermal_state.material + if isinstance(material, Air): + # Air material always uses TPG (via nasa_9_coefficients or thermally_perfect_gas) + uses_tpg = True + + if uses_tpg: + raise ValueError( + "ThermallyPerfectGas model is not supported with the CompressibleIsentropic solver. " + "The CompressibleIsentropic solver uses a 4x4 system that decouples the energy equation " + "and does not support temperature-dependent gas properties. " + "Please use type_name='Compressible' in NavierStokesSolver for thermally perfect gas simulations." + ) + + return params diff --git a/tests/simulation/params/test_validators_material.py b/tests/simulation/params/test_validators_material.py new file mode 100644 index 000000000..44db1adfe --- /dev/null +++ b/tests/simulation/params/test_validators_material.py @@ -0,0 +1,338 @@ +"""Tests for material model validators (NASA9 coefficients, ThermallyPerfectGas).""" + +import pytest + +import flow360 as fl +from flow360.component.simulation.models.material import ( + NASA9CoefficientSet, + NASA9Coefficients, + FrozenSpecies, + ThermallyPerfectGas, +) +from flow360.component.simulation.unit_system import SI_unit_system + + +# ============================================================================= +# NASA9CoefficientSet Tests +# ============================================================================= + + +def test_nasa9_coefficient_set_valid(): + """Test creating a valid NASA9CoefficientSet with exactly 9 coefficients.""" + with SI_unit_system: + coeff_set = 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], + ) + assert len(coeff_set.coefficients) == 9 + + +def test_nasa9_coefficient_set_wrong_count_raises(): + """Test that providing wrong number of coefficients raises ValueError.""" + with SI_unit_system: + with pytest.raises(ValueError, match="requires exactly 9 coefficients"): + NASA9CoefficientSet( + temperature_range_min=200.0 * fl.u.K, + temperature_range_max=1000.0 * fl.u.K, + coefficients=[0.0, 0.0, 3.5], # Only 3 coefficients + ) + + +def test_nasa9_coefficient_set_too_many_coefficients_raises(): + """Test that providing too many coefficients raises ValueError.""" + with SI_unit_system: + with pytest.raises(ValueError, match="requires exactly 9 coefficients"): + NASA9CoefficientSet( + temperature_range_min=200.0 * fl.u.K, + temperature_range_max=1000.0 * fl.u.K, + coefficients=[0.0] * 10, # 10 coefficients + ) + + +# ============================================================================= +# NASA9Coefficients Tests +# ============================================================================= + + +def test_nasa9_coefficients_single_range_valid(): + """Test creating NASA9Coefficients with a single temperature range.""" + with SI_unit_system: + coeffs = NASA9Coefficients( + temperature_ranges=[ + 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], + ) + ] + ) + assert len(coeffs.temperature_ranges) == 1 + + +def test_nasa9_coefficients_multiple_ranges_continuous_valid(): + """Test creating NASA9Coefficients with continuous temperature ranges.""" + with SI_unit_system: + coeffs = NASA9Coefficients( + temperature_ranges=[ + 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], + ), + 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], + ), + ] + ) + assert len(coeffs.temperature_ranges) == 2 + + +def test_nasa9_coefficients_discontinuous_ranges_raises(): + """Test that discontinuous temperature ranges raise ValueError.""" + with SI_unit_system: + with pytest.raises(ValueError, match="Temperature ranges must be continuous"): + NASA9Coefficients( + temperature_ranges=[ + 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], + ), + NASA9CoefficientSet( + temperature_range_min=1500.0 * fl.u.K, # Gap: 1000-1500 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], + ), + ] + ) + + +# ============================================================================= +# ThermallyPerfectGas Tests +# ============================================================================= + + +def _make_species(name: str, mass_fraction: float, t_min: float = 200.0, t_max: float = 6000.0): + """Helper to create a FrozenSpecies with given mass fraction. + + Note: Must be called within a SI_unit_system context. + """ + return FrozenSpecies( + name=name, + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=t_min * fl.u.K, + temperature_range_max=t_max * fl.u.K, + coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ) + ] + ), + mass_fraction=mass_fraction, + ) + + +def test_thermally_perfect_gas_single_species_valid(): + """Test creating ThermallyPerfectGas with a single species (mass_fraction=1.0).""" + with SI_unit_system: + tpg = ThermallyPerfectGas(species=[_make_species("N2", 1.0)]) + assert len(tpg.species) == 1 + assert tpg.species[0].mass_fraction == 1.0 + + +def test_thermally_perfect_gas_multi_species_valid(): + """Test creating ThermallyPerfectGas with multiple species summing to 1.0.""" + with SI_unit_system: + tpg = ThermallyPerfectGas( + species=[ + _make_species("N2", 0.7555), + _make_species("O2", 0.2316), + _make_species("Ar", 0.0129), + ] + ) + assert len(tpg.species) == 3 + total = sum(s.mass_fraction for s in tpg.species) + assert 0.999 <= total <= 1.001 + + +def test_thermally_perfect_gas_mass_fractions_not_sum_to_one_raises(): + """Test that mass fractions not summing to 1.0 raise ValueError.""" + with SI_unit_system: + with pytest.raises(ValueError, match="Mass fractions must sum to 1.0"): + ThermallyPerfectGas( + species=[ + _make_species("N2", 0.5), + _make_species("O2", 0.3), + # Missing 0.2 to sum to 1.0 + ] + ) + + +def test_thermally_perfect_gas_mass_fractions_exceed_one_raises(): + """Test that mass fractions exceeding 1.0 raise ValueError.""" + with SI_unit_system: + with pytest.raises(ValueError, match="Mass fractions must sum to 1.0"): + ThermallyPerfectGas( + species=[ + _make_species("N2", 0.8), + _make_species("O2", 0.4), # Sum = 1.2 + ] + ) + + +def test_thermally_perfect_gas_temperature_ranges_match_valid(): + """Test that species with matching temperature ranges are accepted.""" + with SI_unit_system: + # Both species have same temperature range: 200-1000K, 1000-6000K + n2 = FrozenSpecies( + name="N2", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + 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], + ), + 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], + ), + ] + ), + mass_fraction=0.75, + ) + o2 = FrozenSpecies( + name="O2", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + 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], + ), + 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], + ), + ] + ), + mass_fraction=0.25, + ) + tpg = ThermallyPerfectGas(species=[n2, o2]) + assert len(tpg.species) == 2 + + +def test_thermally_perfect_gas_temperature_ranges_count_mismatch_raises(): + """Test that species with different number of temperature ranges raise ValueError.""" + with SI_unit_system: + # N2 has 2 ranges, O2 has 1 range + n2 = FrozenSpecies( + name="N2", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + 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], + ), + 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], + ), + ] + ), + mass_fraction=0.75, + ) + o2 = FrozenSpecies( + name="O2", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + 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], + ), + ] + ), + mass_fraction=0.25, + ) + with pytest.raises(ValueError, match="same number of temperature ranges"): + ThermallyPerfectGas(species=[n2, o2]) + + +def test_thermally_perfect_gas_temperature_ranges_boundary_mismatch_raises(): + """Test that species with mismatched temperature boundaries raise ValueError.""" + with SI_unit_system: + # N2: 200-1000K, O2: 200-1200K (different boundary) + n2 = FrozenSpecies( + name="N2", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + 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], + ), + ] + ), + mass_fraction=0.75, + ) + o2 = FrozenSpecies( + name="O2", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * fl.u.K, + temperature_range_max=1200.0 * fl.u.K, # Different max + coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ), + ] + ), + mass_fraction=0.25, + ) + with pytest.raises(ValueError, match="boundaries mismatch"): + ThermallyPerfectGas(species=[n2, o2]) + + +# ============================================================================= +# Air with ThermallyPerfectGas Tests +# ============================================================================= + + +def test_air_with_thermally_perfect_gas_valid(): + """Test creating Air material with ThermallyPerfectGas.""" + with SI_unit_system: + tpg = ThermallyPerfectGas( + species=[ + _make_species("N2", 0.7555), + _make_species("O2", 0.2316), + _make_species("Ar", 0.0129), + ] + ) + air = fl.Air(thermally_perfect_gas=tpg) + assert air.thermally_perfect_gas is not None + assert len(air.thermally_perfect_gas.species) == 3 + + +def test_air_with_nasa9_coefficients_valid(): + """Test creating Air material with custom NASA9 coefficients.""" + with SI_unit_system: + air = fl.Air( + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * fl.u.K, + temperature_range_max=6000.0 * fl.u.K, + coefficients=[0.0, 0.0, 3.5, 1e-4, 0.0, 0.0, 0.0, 0.0, 0.0], + ) + ] + ) + ) + assert air.nasa_9_coefficients is not None + assert len(air.nasa_9_coefficients.temperature_ranges) == 1 diff --git a/tests/simulation/translator/test_solver_translator.py b/tests/simulation/translator/test_solver_translator.py index fb03c642e..0c880778d 100644 --- a/tests/simulation/translator/test_solver_translator.py +++ b/tests/simulation/translator/test_solver_translator.py @@ -1467,3 +1467,497 @@ def test_ghost_periodic(): translate_and_compare( processed_params, mesh_unit=1 * u.m, ref_json_file="Flow360_ghost_periodic.json", debug=True ) + + +# ============================================================================= +# NASA 9 Coefficient Translation Tests +# ============================================================================= + + +def test_translate_nasa9_coefficients_single_range(): + """Test translation of NASA9 coefficients with a single temperature range.""" + from flow360.component.simulation.models.material import ( + NASA9CoefficientSet, + NASA9Coefficients, + ) + from flow360.component.simulation.translator.solver_translator import ( + translate_nasa9_coefficients, + ) + + with SI_unit_system: + # Coefficients for constant gamma=1.4 (cp/R = 3.5) + nasa_coeffs = NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ) + ] + ) + + # Translate with reference temperature of 300 K + result = translate_nasa9_coefficients(nasa_coeffs, reference_temperature=300.0) + + assert "temperatureRanges" in result + assert len(result["temperatureRanges"]) == 1 + + range_data = result["temperatureRanges"][0] + # Temperature ranges should be non-dimensionalized: T_nd = T / T_ref + assert range_data["temperatureRangeMin"] == pytest.approx(200.0 / 300.0) + assert range_data["temperatureRangeMax"] == pytest.approx(6000.0 / 300.0) + + # a2 (constant term) should not change + assert range_data["coefficients"][2] == pytest.approx(3.5) + + +def test_translate_nasa9_coefficients_temperature_scaling(): + """Test that temperature-dependent coefficients are scaled correctly.""" + from flow360.component.simulation.models.material import ( + NASA9CoefficientSet, + NASA9Coefficients, + ) + from flow360.component.simulation.translator.solver_translator import ( + translate_nasa9_coefficients, + ) + + T_ref = 300.0 + # cp/R = 3.5 + 0.001*T (temperature dependent) + with SI_unit_system: + nasa_coeffs = NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=[0.0, 0.0, 3.5, 0.001, 0.0, 0.0, 0.0, 0.0, 0.0], + ) + ] + ) + + result = translate_nasa9_coefficients(nasa_coeffs, reference_temperature=T_ref) + coeffs = result["temperatureRanges"][0]["coefficients"] + + # a2 (constant term): no scaling + assert coeffs[2] == pytest.approx(3.5) + + # a3 (T^1 term): scale by T_ref + # a3_nd = a3_dim * T_ref = 0.001 * 300 = 0.3 + assert coeffs[3] == pytest.approx(0.001 * T_ref) + + +def test_translate_thermally_perfect_gas_single_species(): + """Test translation of TPG with a single species (should match direct NASA9 translation).""" + from flow360.component.simulation.models.material import ( + FrozenSpecies, + NASA9CoefficientSet, + NASA9Coefficients, + ThermallyPerfectGas, + ) + from flow360.component.simulation.translator.solver_translator import ( + translate_nasa9_coefficients, + translate_thermally_perfect_gas, + ) + + with SI_unit_system: + coeffs_list = [0.0, 0.0, 3.5, 0.001, 0.0, 0.0, 0.0, 0.0, 0.0] + nasa_coeffs = NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=coeffs_list, + ) + ] + ) + tpg = ThermallyPerfectGas( + species=[ + FrozenSpecies( + name="Air", + nasa_9_coefficients=nasa_coeffs, + mass_fraction=1.0, + ) + ] + ) + + T_ref = 300.0 + result_tpg = translate_thermally_perfect_gas(tpg, reference_temperature=T_ref) + result_single = translate_nasa9_coefficients(nasa_coeffs, reference_temperature=T_ref) + + # Single species with mass_fraction=1.0 should give same result as direct translation + assert len(result_tpg["temperatureRanges"]) == len(result_single["temperatureRanges"]) + + for i in range(len(result_tpg["temperatureRanges"])): + tpg_coeffs = result_tpg["temperatureRanges"][i]["coefficients"] + single_coeffs = result_single["temperatureRanges"][i]["coefficients"] + for j in range(9): + assert tpg_coeffs[j] == pytest.approx(single_coeffs[j]) + + +def test_translate_thermally_perfect_gas_mass_fraction_weighting(): + """Test that multi-species coefficients are correctly mass-fraction weighted.""" + from flow360.component.simulation.models.material import ( + FrozenSpecies, + NASA9CoefficientSet, + NASA9Coefficients, + ThermallyPerfectGas, + ) + from flow360.component.simulation.translator.solver_translator import ( + translate_thermally_perfect_gas, + ) + + # Species A: cp/R = 3.0 (a2 = 3.0) + # Species B: cp/R = 4.0 (a2 = 4.0) + # Mass fractions: A = 0.6, B = 0.4 + # Expected combined: cp/R = 0.6 * 3.0 + 0.4 * 4.0 = 1.8 + 1.6 = 3.4 + + with SI_unit_system: + species_a = FrozenSpecies( + name="A", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=[0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ) + ] + ), + mass_fraction=0.6, + ) + species_b = FrozenSpecies( + name="B", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=[0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ) + ] + ), + mass_fraction=0.4, + ) + tpg = ThermallyPerfectGas(species=[species_a, species_b]) + + result = translate_thermally_perfect_gas(tpg, reference_temperature=300.0) + coeffs = result["temperatureRanges"][0]["coefficients"] + + # a2 should be mass-fraction weighted: 0.6 * 3.0 + 0.4 * 4.0 = 3.4 + assert coeffs[2] == pytest.approx(3.4) + + +def test_translate_thermally_perfect_gas_all_coefficients_weighted(): + """Test that all 9 coefficients are correctly mass-fraction weighted.""" + from flow360.component.simulation.models.material import ( + FrozenSpecies, + NASA9CoefficientSet, + NASA9Coefficients, + ThermallyPerfectGas, + ) + from flow360.component.simulation.translator.solver_translator import ( + translate_thermally_perfect_gas, + ) + + # Two species with different coefficients + coeffs_a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0] + coeffs_b = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0] + y_a, y_b = 0.7, 0.3 + + with SI_unit_system: + species_a = FrozenSpecies( + name="A", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=coeffs_a, + ) + ] + ), + mass_fraction=y_a, + ) + species_b = FrozenSpecies( + name="B", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=coeffs_b, + ) + ] + ), + mass_fraction=y_b, + ) + tpg = ThermallyPerfectGas(species=[species_a, species_b]) + + T_ref = 300.0 + result = translate_thermally_perfect_gas(tpg, reference_temperature=T_ref) + result_coeffs = result["temperatureRanges"][0]["coefficients"] + + # Compute expected combined coefficients (before non-dimensionalization) + expected_combined = [y_a * coeffs_a[i] + y_b * coeffs_b[i] for i in range(9)] + + # Apply non-dimensionalization scaling factors + t_scale = T_ref + t_scale_inv = 1.0 / T_ref + expected_nd = [ + expected_combined[0] * t_scale_inv**2, # a0 (T^-2) + expected_combined[1] * t_scale_inv, # a1 (T^-1) + expected_combined[2], # a2 (constant) + expected_combined[3] * t_scale, # a3 (T^1) + expected_combined[4] * t_scale**2, # a4 (T^2) + expected_combined[5] * t_scale**3, # a5 (T^3) + expected_combined[6] * t_scale**4, # a6 (T^4) + expected_combined[7] * t_scale_inv, # a7 (enthalpy, 1/T) + expected_combined[8], # a8 (entropy, constant) + ] + + for i in range(9): + assert result_coeffs[i] == pytest.approx(expected_nd[i]), f"Coefficient a{i} mismatch" + + +def test_translate_thermally_perfect_gas_multiple_temperature_ranges(): + """Test translation with multiple temperature ranges.""" + from flow360.component.simulation.models.material import ( + FrozenSpecies, + NASA9CoefficientSet, + NASA9Coefficients, + ThermallyPerfectGas, + ) + from flow360.component.simulation.translator.solver_translator import ( + translate_thermally_perfect_gas, + ) + + # Both species have two temperature ranges with matching boundaries + with SI_unit_system: + species_a = FrozenSpecies( + name="A", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=1000.0 * u.K, + coefficients=[0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ), + NASA9CoefficientSet( + temperature_range_min=1000.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=[0.0, 0.0, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ), + ] + ), + mass_fraction=0.75, + ) + species_b = FrozenSpecies( + name="B", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=1000.0 * u.K, + coefficients=[0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ), + NASA9CoefficientSet( + temperature_range_min=1000.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=[0.0, 0.0, 4.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ), + ] + ), + mass_fraction=0.25, + ) + tpg = ThermallyPerfectGas(species=[species_a, species_b]) + + result = translate_thermally_perfect_gas(tpg, reference_temperature=300.0) + + assert len(result["temperatureRanges"]) == 2 + + # First range: 0.75 * 3.0 + 0.25 * 4.0 = 2.25 + 1.0 = 3.25 + assert result["temperatureRanges"][0]["coefficients"][2] == pytest.approx(3.25) + + # Second range: 0.75 * 3.5 + 0.25 * 4.5 = 2.625 + 1.125 = 3.75 + assert result["temperatureRanges"][1]["coefficients"][2] == pytest.approx(3.75) + + +def test_nasa9_a7_correction_ensures_correct_internal_energy(): + """ + Test that the a7 coefficient correction ensures correct internal energy at T_nd=1. + + The NASA database a7 coefficient is calibrated for absolute dimensional enthalpy, + but Flow360's non-dimensionalization expects e = cv*T at T_ref. + + The correction computes gamma from the polynomial itself (not a fixed value), + ensuring internal energy consistency: e/R = cv/R * T = (1/(gamma-1)) * T at T_nd=1. + """ + import math + + from flow360.component.simulation.models.material import ( + NASA9CoefficientSet, + NASA9Coefficients, + ) + from flow360.component.simulation.translator.solver_translator import ( + translate_nasa9_coefficients, + ) + + T_ref = 300.0 # K + + # Use real NASA-like coefficients with a non-zero a7 that would cause issues + # This is similar to the combined air coefficients used in the ShockTube test + dim_coeffs = [ + 8765.75, # a0 + -176.23, # a1 + 4.887, # a2 (cp/R ~ 3.5 for diatomic gas at low T) + -0.00545, # a3 + 1.03e-05, # a4 + -7.74e-09, # a5 + 2.14e-12, # a6 + -1585.05, # a7 - Original NASA value that causes problems + -3.57, # a8 + ] + + nasa_coeffs = NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=dim_coeffs, + ) + ] + ) + + # Use SI_unit_system context for the test + with SI_unit_system: + result = translate_nasa9_coefficients(nasa_coeffs, T_ref) + + # Get the translated (non-dimensionalized) coefficients + nd_coeffs = result["temperatureRanges"][0]["coefficients"] + + # Compute h/R at T_nd=1 using the translated polynomial + T_nd = 1.0 + h_over_R = ( + -nd_coeffs[0] / T_nd + + nd_coeffs[1] * math.log(T_nd) + + nd_coeffs[2] * T_nd + + (nd_coeffs[3] / 2) * T_nd**2 + + (nd_coeffs[4] / 3) * T_nd**3 + + (nd_coeffs[5] / 4) * T_nd**4 + + (nd_coeffs[6] / 5) * T_nd**5 + + nd_coeffs[7] + ) + + # Compute e/R = h/R - T (non-dimensional) + e_over_R = h_over_R - T_nd + + # Compute expected e/R from the polynomial's own gamma at T_ref + # cp/R at T_ref (dimensional) + cp_over_R = ( + dim_coeffs[0] * T_ref ** (-2) + + dim_coeffs[1] * T_ref ** (-1) + + dim_coeffs[2] + + dim_coeffs[3] * T_ref + + dim_coeffs[4] * T_ref**2 + + dim_coeffs[5] * T_ref**3 + + dim_coeffs[6] * T_ref**4 + ) + # cv/R = cp/R - 1 + cv_over_R = cp_over_R - 1 + # gamma = cp/cv + gamma = cp_over_R / cv_over_R + # Expected e/R at T_nd=1 is cv/R * T = (1/(gamma-1)) * 1 + expected_e_over_R = 1.0 / (gamma - 1) + + assert e_over_R == pytest.approx( + expected_e_over_R, rel=1e-6 + ), f"Internal energy e/R={e_over_R} at T_nd=1, expected {expected_e_over_R} (gamma={gamma})" + + +def test_nasa9_a7_correction_with_thermally_perfect_gas(): + """ + Test the a7 correction through the translate_thermally_perfect_gas function. + + This tests the full translation pipeline with multi-species TPG. + """ + import math + + from flow360.component.simulation.models.material import ( + FrozenSpecies, + NASA9CoefficientSet, + NASA9Coefficients, + ThermallyPerfectGas, + ) + from flow360.component.simulation.translator.solver_translator import ( + translate_thermally_perfect_gas, + ) + + T_ref = 300.0 # K + + # Use real NASA-like coefficients + dim_coeffs = [ + 8765.75, # a0 + -176.23, # a1 + 4.887, # a2 + -0.00545, # a3 + 1.03e-05, # a4 + -7.74e-09, # a5 + 2.14e-12, # a6 + -1585.05, # a7 - Original NASA value + -3.57, # a8 + ] + + # Create a simple single-species TPG with real-like NASA coefficients + with SI_unit_system: + species = FrozenSpecies( + name="air", + nasa_9_coefficients=NASA9Coefficients( + temperature_ranges=[ + NASA9CoefficientSet( + temperature_range_min=200.0 * u.K, + temperature_range_max=6000.0 * u.K, + coefficients=dim_coeffs, + ) + ] + ), + mass_fraction=1.0, + ) + tpg = ThermallyPerfectGas(species=[species]) + + result = translate_thermally_perfect_gas(tpg, reference_temperature=T_ref) + + # Get the translated coefficients + nd_coeffs = result["temperatureRanges"][0]["coefficients"] + + # Compute h/R at T_nd=1 + T_nd = 1.0 + h_over_R = ( + -nd_coeffs[0] / T_nd + + nd_coeffs[1] * math.log(T_nd) + + nd_coeffs[2] * T_nd + + (nd_coeffs[3] / 2) * T_nd**2 + + (nd_coeffs[4] / 3) * T_nd**3 + + (nd_coeffs[5] / 4) * T_nd**4 + + (nd_coeffs[6] / 5) * T_nd**5 + + nd_coeffs[7] + ) + + # e/R = h/R - T + e_over_R = h_over_R - T_nd + + # Compute expected e/R from the polynomial's own gamma at T_ref + cp_over_R = ( + dim_coeffs[0] * T_ref ** (-2) + + dim_coeffs[1] * T_ref ** (-1) + + dim_coeffs[2] + + dim_coeffs[3] * T_ref + + dim_coeffs[4] * T_ref**2 + + dim_coeffs[5] * T_ref**3 + + dim_coeffs[6] * T_ref**4 + ) + cv_over_R = cp_over_R - 1 + gamma = cp_over_R / cv_over_R + expected_e_over_R = 1.0 / (gamma - 1) + + assert e_over_R == pytest.approx( + expected_e_over_R, rel=1e-6 + ), f"Internal energy e/R={e_over_R} at T_nd=1, expected {expected_e_over_R} (gamma={gamma})"