diff --git a/flow360/__init__.py b/flow360/__init__.py index dea86b3e3..f397a3322 100644 --- a/flow360/__init__.py +++ b/flow360/__init__.py @@ -50,6 +50,7 @@ SpalartAllmaras, SpalartAllmarasModelConstants, TransitionModelSolver, + TurbulenceModelControls, ) from flow360.component.simulation.models.surface_models import ( Freestream, @@ -248,6 +249,7 @@ "Box", "GenericReferenceCondition", "TransitionModelSolver", + "TurbulenceModelControls", "Pressure", "TotalPressure", "Rotational", diff --git a/flow360/component/simulation/models/solver_numerics.py b/flow360/component/simulation/models/solver_numerics.py index f5b9f0fbd..cb3f2c001 100644 --- a/flow360/component/simulation/models/solver_numerics.py +++ b/flow360/component/simulation/models/solver_numerics.py @@ -10,7 +10,7 @@ from __future__ import annotations from abc import ABCMeta -from typing import Annotated, Dict, Literal, Optional, Union +from typing import Annotated, Dict, List, Literal, Optional, Union import numpy as np import pydantic as pd @@ -22,7 +22,7 @@ Flow360BaseModel, ) from flow360.component.simulation.framework.entity_base import EntityList -from flow360.component.simulation.primitives import Box +from flow360.component.simulation.primitives import Box, GenericVolume # from .time_stepping import UnsteadyTimeStepping @@ -202,6 +202,48 @@ class KOmegaSSTModelConstants(Flow360BaseModel): ] +class TurbulenceModelControls(Flow360BaseModel): + """ + :class:`TurbulenceModelControls` class specifies modeling constants and enforces turbulence model + behavior on a zonal basis, as defined by mesh entities or boxes in space. These controls + supersede the global turbulence model solver settings. + + Example + _______ + >>> fl.TurbulenceModelControls( + ... modeling_constants=fl.SpalartAllmarasConstants(C_w2=2.718), + ... enforcement="RANS", + ... entities=[ + ... volume_mesh["block-1"], + ... fl.Box.from_principal_axes( + ... name="box", + ... axes=[(0, 1, 0), (0, 0, 1)], + ... center=(0, 0, 0) * fl.u.m, + ... size=(0.2, 0.3, 2) * fl.u.m, + ... ), + ... ], + ... ) + """ + + modeling_constants: Optional[TurbulenceModelConstants] = pd.Field( + None, + description="A class of `SpalartAllmarasModelConstants` or `KOmegaSSTModelConstants` used to " + + "specify constants in specific regions of the domain.", + ) + + enforcement: Optional[Literal["RANS", "LES"]] = pd.Field( + None, description="Force RANS or LES mode in a specific control region." + ) + + entities: EntityList[GenericVolume, Box] = pd.Field( + alias="volumes", + description="The entity in which to apply the `TurbulenceMOdelControls``. " + + "The entity should be defined by :class:`Box` or zones from the geometry/volume mesh." + + "The axes of entity must be specified to serve as the the principle axes of the " + + "`TurbulenceModelControls` region.", + ) + + class DetachedEddySimulation(Flow360BaseModel): """ :class:`DetachedEddySimulation` class is used for running hybrid RANS-LES simulations @@ -295,6 +337,38 @@ class TurbulenceModelSolver(GenericSolverSettings, metaclass=ABCMeta): False, description="Rotation correction for the turbulence model." ) + controls: Optional[List[TurbulenceModelControls]] = pd.Field( + None, + strict=True, # Note: To ensure propoer err msg when none-list is fed. + description="List of control zones to enforce specific turbulence model constants " + + "and behavior.", + ) + + @pd.model_validator(mode="after") + def _check_zonal_modeling_constants_consistency(self) -> Self: + if self.controls is None: + return self + + for index, control in enumerate(self.controls): + if control.modeling_constants is None: + continue + if not isinstance( + control.modeling_constants, SpalartAllmarasModelConstants + ) and isinstance(self, SpalartAllmaras): + raise ValueError( + "Turbulence model is SpalartAllmaras, but controls.modeling" + "_constants is of a conflicting class " + f"in control region {index}." + ) + if not isinstance(control.modeling_constants, KOmegaSSTModelConstants) and isinstance( + self, KOmegaSST + ): + raise ValueError( + "Turbulence model is KOmegaSST, but controls.modeling_constants" + f" is of a conflicting class in control region {index}." + ) + return self + class KOmegaSST(TurbulenceModelSolver): """ diff --git a/flow360/component/simulation/simulation_params.py b/flow360/component/simulation/simulation_params.py index fda8967aa..ae049056e 100644 --- a/flow360/component/simulation/simulation_params.py +++ b/flow360/component/simulation/simulation_params.py @@ -77,6 +77,7 @@ _check_consistency_wall_function_and_surface_output, _check_duplicate_entities_in_models, _check_duplicate_isosurface_names, + _check_hybrid_model_to_use_zonal_enforcement, _check_low_mach_preconditioner_output, _check_numerical_dissipation_factor_output, _check_parent_volume_is_rotating, @@ -446,6 +447,11 @@ def check_unsteadiness_to_use_hybrid_model(self): """Only allow hybrid RANS-LES output field for unsteady simulations""" return _check_unsteadiness_to_use_hybrid_model(self) + @pd.model_validator(mode="after") + def check_hybrid_model_to_use_zonal_enforcement(self): + """Only allow LES/RANS zonal enforcement in hybrid RANS-LES mode""" + return _check_hybrid_model_to_use_zonal_enforcement(self) + @pd.model_validator(mode="after") def check_unsteadiness_to_use_aero_acoustics(self): """Only allow Aero acoustics when using unsteady simulation""" diff --git a/flow360/component/simulation/translator/solver_translator.py b/flow360/component/simulation/translator/solver_translator.py index 33d450e89..7ab6c2402 100644 --- a/flow360/component/simulation/translator/solver_translator.py +++ b/flow360/component/simulation/translator/solver_translator.py @@ -961,6 +961,49 @@ def get_navier_stokes_initial_condition( return initial_condition_dict +def rename_modeling_constants(modeling_constants): + """Rename the modeling constants to what the solver reads""" + if modeling_constants.get("typeName", None) == "SpalartAllmarasConsts": + replace_dict_key(modeling_constants, "CDES", "C_DES") + replace_dict_key(modeling_constants, "CD", "C_d") + replace_dict_key(modeling_constants, "CCb1", "C_cb1") + replace_dict_key(modeling_constants, "CCb2", "C_cb2") + replace_dict_key(modeling_constants, "CSigma", "C_sigma") + replace_dict_key(modeling_constants, "CV1", "C_v1") + replace_dict_key(modeling_constants, "CVonKarman", "C_vonKarman") + replace_dict_key(modeling_constants, "CW2", "C_w2") + replace_dict_key(modeling_constants, "CT3", "C_t3") + replace_dict_key(modeling_constants, "CT4", "C_t4") + replace_dict_key(modeling_constants, "CMinRd", "C_min_rd") + + if modeling_constants.get("typeName", None) == "kOmegaSSTConsts": + replace_dict_key(modeling_constants, "CDES1", "C_DES1") + replace_dict_key(modeling_constants, "CDES2", "C_DES2") + replace_dict_key(modeling_constants, "CD1", "C_d1") + replace_dict_key(modeling_constants, "CD2", "C_d2") + replace_dict_key(modeling_constants, "CSigmaK1", "C_sigma_k1") + replace_dict_key(modeling_constants, "CSigmaK2", "C_sigma_k2") + replace_dict_key(modeling_constants, "CSigmaOmega1", "C_sigma_omega1") + replace_dict_key(modeling_constants, "CSigmaOmega2", "C_sigma_omega2") + replace_dict_key(modeling_constants, "CAlpha1", "C_alpha1") + replace_dict_key(modeling_constants, "CBeta1", "C_beta1") + replace_dict_key(modeling_constants, "CBeta2", "C_beta2") + replace_dict_key(modeling_constants, "CBetaStar", "C_beta_star") + + modeling_constants.pop("typeName") # Not read by solver + + +def update_controls_modeling_constants(controls, translated): + """Upading the modelingConstants entries for each control""" + if controls is not None: + for control in translated["turbulenceModelSolver"]["controls"]: + control_modeling_constants = control.get("modelingConstants", None) + if control_modeling_constants is None: + continue + rename_modeling_constants(control_modeling_constants) + control["modelConstants"] = control.pop("modelingConstants") + + # pylint: disable=too-many-statements # pylint: disable=too-many-branches # pylint: disable=too-many-locals @@ -1074,34 +1117,7 @@ def get_solver_json( replace_dict_key(translated["turbulenceModelSolver"], "typeName", "modelType") modeling_constants = translated["turbulenceModelSolver"].get("modelingConstants", None) if modeling_constants is not None: - if modeling_constants.get("typeName", None) == "SpalartAllmarasConsts": - replace_dict_key(modeling_constants, "CDES", "C_DES") - replace_dict_key(modeling_constants, "CD", "C_d") - replace_dict_key(modeling_constants, "CCb1", "C_cb1") - replace_dict_key(modeling_constants, "CCb2", "C_cb2") - replace_dict_key(modeling_constants, "CSigma", "C_sigma") - replace_dict_key(modeling_constants, "CV1", "C_v1") - replace_dict_key(modeling_constants, "CVonKarman", "C_vonKarman") - replace_dict_key(modeling_constants, "CW2", "C_w2") - replace_dict_key(modeling_constants, "CT3", "C_t3") - replace_dict_key(modeling_constants, "CT4", "C_t4") - replace_dict_key(modeling_constants, "CMinRd", "C_min_rd") - - if modeling_constants.get("typeName", None) == "kOmegaSSTConsts": - replace_dict_key(modeling_constants, "CDES1", "C_DES1") - replace_dict_key(modeling_constants, "CDES2", "C_DES2") - replace_dict_key(modeling_constants, "CD1", "C_d1") - replace_dict_key(modeling_constants, "CD2", "C_d2") - replace_dict_key(modeling_constants, "CSigmaK1", "C_sigma_k1") - replace_dict_key(modeling_constants, "CSigmaK2", "C_sigma_k2") - replace_dict_key(modeling_constants, "CSigmaOmega1", "C_sigma_omega1") - replace_dict_key(modeling_constants, "CSigmaOmega2", "C_sigma_omega2") - replace_dict_key(modeling_constants, "CAlpha1", "C_alpha1") - replace_dict_key(modeling_constants, "CBeta1", "C_beta1") - replace_dict_key(modeling_constants, "CBeta2", "C_beta2") - replace_dict_key(modeling_constants, "CBetaStar", "C_beta_star") - - modeling_constants.pop("typeName") # Not read by solver + rename_modeling_constants(modeling_constants) translated["turbulenceModelSolver"]["modelConstants"] = translated[ "turbulenceModelSolver" ].pop("modelingConstants") @@ -1124,6 +1140,10 @@ def get_solver_json( translated["turbulenceModelSolver"]["ZDES"] = False translated["turbulenceModelSolver"]["gridSizeForLES"] = "maxEdgeLength" + update_controls_modeling_constants( + model.turbulence_model_solver.controls, translated + ) + if not isinstance(model.transition_model_solver, NoneSolver): # baseline dictionary dump for transition model object translated["transitionModelSolver"] = dump_dict(model.transition_model_solver) diff --git a/flow360/component/simulation/validation/validation_simulation_params.py b/flow360/component/simulation/validation/validation_simulation_params.py index e9b0b8c60..6f9b5179d 100644 --- a/flow360/component/simulation/validation/validation_simulation_params.py +++ b/flow360/component/simulation/validation/validation_simulation_params.py @@ -225,6 +225,30 @@ def _check_unsteadiness_to_use_hybrid_model(v): return v +def _check_hybrid_model_to_use_zonal_enforcement(v): + models = v.models + if not models: + return v + + for model in models: + if isinstance(model, Fluid): + turbulence_model_solver = model.turbulence_model_solver + if not isinstance(turbulence_model_solver, NoneSolver): + if turbulence_model_solver.controls is None: + continue + for index, control in enumerate(turbulence_model_solver.controls): + if ( + control.enforcement is not None + and turbulence_model_solver.hybrid_model is None + ): + raise ValueError( + f"Control region {index} must be running in hybrid RANS-LES mode to " + "apply zonal turbulence enforcement." + ) + + return v + + def _check_cht_solver_settings(params): has_heat_transfer = False diff --git a/tests/simulation/converter/ref/ref_monitor.json b/tests/simulation/converter/ref/ref_monitor.json index 834d23adb..bb2931d02 100644 --- a/tests/simulation/converter/ref/ref_monitor.json +++ b/tests/simulation/converter/ref/ref_monitor.json @@ -1 +1 @@ -{"version":"25.6.0b1","unit_system":{"name":"SI"},"meshing":null,"reference_geometry":null,"operating_condition":null,"models":[{"material":{"type":"air","name":"air","dynamic_viscosity":{"reference_viscosity":{"value":0.00001716,"units":"Pa*s"},"reference_temperature":{"value":273.15,"units":"K"},"effective_temperature":{"value":110.4,"units":"K"}}},"initial_condition":{"type_name":"NavierStokesInitialCondition","constants":null,"rho":"rho","u":"u","v":"v","w":"w","p":"p"},"type":"Fluid","navier_stokes_solver":{"absolute_tolerance":1e-10,"relative_tolerance":0.0,"order_of_accuracy":2,"equation_evaluation_frequency":1,"linear_solver":{"max_iterations":30,"absolute_tolerance":null,"relative_tolerance":null},"private_attribute_dict":null,"CFL_multiplier":1.0,"kappa_MUSCL":-1.0,"numerical_dissipation_factor":1.0,"limit_velocity":false,"limit_pressure_density":false,"type_name":"Compressible","low_mach_preconditioner":false,"low_mach_preconditioner_threshold":null,"update_jacobian_frequency":4,"max_force_jac_update_physical_steps":0},"turbulence_model_solver":{"absolute_tolerance":1e-8,"relative_tolerance":0.0,"order_of_accuracy":2,"equation_evaluation_frequency":4,"linear_solver":{"max_iterations":20,"absolute_tolerance":null,"relative_tolerance":null},"private_attribute_dict":null,"CFL_multiplier":2.0,"type_name":"SpalartAllmaras","reconstruction_gradient_limiter":0.5,"quadratic_constitutive_relation":false,"modeling_constants":{"type_name":"SpalartAllmarasConsts","C_DES":0.72,"C_d":8.0,"C_cb1":0.1355,"C_cb2":0.622,"C_sigma":0.6666666666666666,"C_v1":7.1,"C_vonKarman":0.41,"C_w2":0.3,"C_t3":1.2,"C_t4":0.5,"C_min_rd":10.0},"update_jacobian_frequency":4,"max_force_jac_update_physical_steps":0,"hybrid_model":null,"rotation_correction":false},"transition_model_solver":{"type_name":"None"}}],"time_stepping":{"type_name":"Steady","max_steps":2000,"CFL":{"type":"adaptive","min":0.1,"max":10000.0,"max_relative_change":1.0,"convergence_limiting_factor":0.25}},"user_defined_dynamics":null,"user_defined_fields":[],"outputs":[{"name":"R1","entities":{"stored_entities":[{"private_attribute_registry_bucket_name":"PointEntityType","private_attribute_entity_type_name":"Point","private_attribute_id":"b9de2bce-36c1-4bbf-af0a-2c6a2ab713a4","name":"Point-0","location":{"value":[2.694298,0.0,1.0195910000000001],"units":"m"}}]},"output_fields":{"items":["primitiveVars"]},"output_type":"ProbeOutput"},{"name":"V3","entities":{"stored_entities":[{"private_attribute_registry_bucket_name":"PointEntityType","private_attribute_entity_type_name":"Point","private_attribute_id":"a79cffc0-31d0-499d-906c-f271c2320166","name":"Point-1","location":{"value":[4.007,0.0,-0.31760000000000005],"units":"m"}},{"private_attribute_registry_bucket_name":"PointEntityType","private_attribute_entity_type_name":"Point","private_attribute_id":"8947eb10-fc59-4102-b9c7-168a91ca22b9","name":"Point-2","location":{"value":[4.007,0.0,-0.29760000000000003],"units":"m"}},{"private_attribute_registry_bucket_name":"PointEntityType","private_attribute_entity_type_name":"Point","private_attribute_id":"27ac4e03-592b-4dba-8fa1-8f6678087a96","name":"Point-3","location":{"value":[4.007,0.0,-0.2776],"units":"m"}}]},"output_fields":{"items":["mut"]},"output_type":"ProbeOutput"}],"private_attribute_asset_cache":{"project_length_unit":null,"project_entity_info":null, "use_inhouse_mesher": false, "use_geometry_AI": false}} +{"version":"25.6.0b1","unit_system":{"name":"SI"},"meshing":null,"reference_geometry":null,"operating_condition":null,"models":[{"material":{"type":"air","name":"air","dynamic_viscosity":{"reference_viscosity":{"value":0.00001716,"units":"Pa*s"},"reference_temperature":{"value":273.15,"units":"K"},"effective_temperature":{"value":110.4,"units":"K"}}},"initial_condition":{"type_name":"NavierStokesInitialCondition","constants":null,"rho":"rho","u":"u","v":"v","w":"w","p":"p"},"type":"Fluid","navier_stokes_solver":{"absolute_tolerance":1e-10,"relative_tolerance":0.0,"order_of_accuracy":2,"equation_evaluation_frequency":1,"linear_solver":{"max_iterations":30,"absolute_tolerance":null,"relative_tolerance":null},"private_attribute_dict":null,"CFL_multiplier":1.0,"kappa_MUSCL":-1.0,"numerical_dissipation_factor":1.0,"limit_velocity":false,"limit_pressure_density":false,"type_name":"Compressible","low_mach_preconditioner":false,"low_mach_preconditioner_threshold":null,"update_jacobian_frequency":4,"max_force_jac_update_physical_steps":0},"turbulence_model_solver":{"absolute_tolerance":1e-8,"relative_tolerance":0.0,"order_of_accuracy":2,"equation_evaluation_frequency":4,"linear_solver":{"max_iterations":20,"absolute_tolerance":null,"relative_tolerance":null},"private_attribute_dict":null,"CFL_multiplier":2.0,"type_name":"SpalartAllmaras","reconstruction_gradient_limiter":0.5,"quadratic_constitutive_relation":false,"modeling_constants":{"type_name":"SpalartAllmarasConsts","C_DES":0.72,"C_d":8.0,"C_cb1":0.1355,"C_cb2":0.622,"C_sigma":0.6666666666666666,"C_v1":7.1,"C_vonKarman":0.41,"C_w2":0.3,"C_t3":1.2,"C_t4":0.5,"C_min_rd":10.0},"update_jacobian_frequency":4,"max_force_jac_update_physical_steps":0,"hybrid_model":null,"rotation_correction":false, "controls":null},"transition_model_solver":{"type_name":"None"}}],"time_stepping":{"type_name":"Steady","max_steps":2000,"CFL":{"type":"adaptive","min":0.1,"max":10000.0,"max_relative_change":1.0,"convergence_limiting_factor":0.25}},"user_defined_dynamics":null,"user_defined_fields":[],"outputs":[{"name":"R1","entities":{"stored_entities":[{"private_attribute_registry_bucket_name":"PointEntityType","private_attribute_entity_type_name":"Point","private_attribute_id":"b9de2bce-36c1-4bbf-af0a-2c6a2ab713a4","name":"Point-0","location":{"value":[2.694298,0.0,1.0195910000000001],"units":"m"}}]},"output_fields":{"items":["primitiveVars"]},"output_type":"ProbeOutput"},{"name":"V3","entities":{"stored_entities":[{"private_attribute_registry_bucket_name":"PointEntityType","private_attribute_entity_type_name":"Point","private_attribute_id":"a79cffc0-31d0-499d-906c-f271c2320166","name":"Point-1","location":{"value":[4.007,0.0,-0.31760000000000005],"units":"m"}},{"private_attribute_registry_bucket_name":"PointEntityType","private_attribute_entity_type_name":"Point","private_attribute_id":"8947eb10-fc59-4102-b9c7-168a91ca22b9","name":"Point-2","location":{"value":[4.007,0.0,-0.29760000000000003],"units":"m"}},{"private_attribute_registry_bucket_name":"PointEntityType","private_attribute_entity_type_name":"Point","private_attribute_id":"27ac4e03-592b-4dba-8fa1-8f6678087a96","name":"Point-3","location":{"value":[4.007,0.0,-0.2776],"units":"m"}}]},"output_fields":{"items":["mut"]},"output_type":"ProbeOutput"}],"private_attribute_asset_cache":{"project_length_unit":null,"project_entity_info":null, "use_inhouse_mesher": false, "use_geometry_AI": false}} diff --git a/tests/simulation/params/test_validators_params.py b/tests/simulation/params/test_validators_params.py index 59ac9651e..9bcf2f562 100644 --- a/tests/simulation/params/test_validators_params.py +++ b/tests/simulation/params/test_validators_params.py @@ -20,7 +20,12 @@ from flow360.component.simulation.models.material import SolidMaterial, aluminum from flow360.component.simulation.models.solver_numerics import ( DetachedEddySimulation, + KOmegaSST, + KOmegaSSTModelConstants, + SpalartAllmaras, + SpalartAllmarasModelConstants, TransitionModelSolver, + TurbulenceModelControls, ) from flow360.component.simulation.models.surface_models import ( Freestream, @@ -318,6 +323,78 @@ def test_hybrid_model_for_unsteady_validator( SimulationParams(models=[fluid_model_with_hybrid_model]) +def test_hybrid_model_to_use_zonal_enforcement(fluid_model, fluid_model_with_hybrid_model): + + fluid_model_with_hybrid_model.turbulence_model_solver.controls = [ + TurbulenceModelControls(enforcement="RANS", entities=[GenericVolume(name="block-1")]) + ] + + # Valid simulation params + with SI_unit_system: + params = SimulationParams( + models=[fluid_model_with_hybrid_model], + time_stepping=Unsteady(steps=12, step_size=0.1 * u.s), + ) + + assert params + + fluid_model.turbulence_model_solver.controls = [ + TurbulenceModelControls(enforcement="RANS", entities=[GenericVolume(name="block-1")]) + ] + + message = "Control region 0 must be running in hybrid RANS-LES mode to apply zonal turbulence enforcement." + with SI_unit_system, pytest.raises(ValueError, match=re.escape(message)): + SimulationParams( + models=[fluid_model], + time_stepping=Unsteady(steps=12, step_size=0.1 * u.s), + ) + + +def test_zonal_modeling_constants_consistency(fluid_model_with_hybrid_model): + fluid_model_with_hybrid_model.turbulence_model_solver.controls = [ + TurbulenceModelControls( + enforcement="RANS", + modeling_constants=SpalartAllmarasModelConstants(), + entities=[GenericVolume(name="block-1")], + ) + ] + + # Valid simulation params + with SI_unit_system: + params = SimulationParams( + models=[fluid_model_with_hybrid_model], + time_stepping=Unsteady(steps=12, step_size=0.1 * u.s), + ) + + assert params + + message = "Turbulence model is SpalartAllmaras, but controls.modeling_constants is of a " + "conflicting class in control region 0." + with SI_unit_system, pytest.raises(ValueError, match=re.escape(message)): + TurbulenceModelSolver = SpalartAllmaras( + controls=[ + TurbulenceModelControls( + enforcement="LES", + modeling_constants=KOmegaSSTModelConstants(), + entities=[GenericVolume(name="block-1")], + ) + ] + ) + + message = "Turbulence model is KOmegaSST, but controls.modeling_constants is of a " + "conflicting class in control region 0." + with SI_unit_system, pytest.raises(ValueError, match=re.escape(message)): + TurbulenceModelSolver = KOmegaSST( + controls=[ + TurbulenceModelControls( + enforcement="LES", + modeling_constants=SpalartAllmarasModelConstants(), + entities=[GenericVolume(name="block-1")], + ) + ] + ) + + def test_cht_solver_settings_validator( fluid_model, ):