diff --git a/flow360/component/simulation/outputs/outputs.py b/flow360/component/simulation/outputs/outputs.py
index 558f30daa..76e0369fa 100644
--- a/flow360/component/simulation/outputs/outputs.py
+++ b/flow360/component/simulation/outputs/outputs.py
@@ -37,6 +37,7 @@
     Surface,
 )
 from flow360.component.simulation.unit_system import LengthType
+from flow360.component.simulation.user_code import ValueOrExpression
 from flow360.component.simulation.validation.validation_context import (
     ALL,
     CASE,
@@ -260,7 +261,9 @@ class VolumeOutput(_AnimationAndFileFormatSettings):
     """
 
     name: Optional[str] = pd.Field("Volume output", description="Name of the `VolumeOutput`.")
-    output_fields: UniqueItemList[Union[VolumeFieldNames, str]] = pd.Field(
+    # TODO: Not all SolverVariables can be used here.
+    # TODO: `Expression` and `SolverVariable` still need business logic for validation (Surface Field/Expression?)
+    output_fields: UniqueItemList[Union[VolumeFieldNames, str, ValueOrExpression]] = pd.Field(
         description="List of output variables. Including :ref:`universal output variables<UniversalVariablesV2>`,"
         " :ref:`variables specific to VolumeOutput<VolumeAndSliceSpecificVariablesV2>`"
         " and :class:`UserDefinedField`."
diff --git a/flow360/component/simulation/services.py b/flow360/component/simulation/services.py
index 4ca5d7550..e4ddbf694 100644
--- a/flow360/component/simulation/services.py
+++ b/flow360/component/simulation/services.py
@@ -522,7 +522,7 @@ def _translate_simulation_json(
     translation_func=None,
 ):
     """
-    Get JSON for surface meshing from a given simulaiton JSON.
+    Get JSON for surface meshing from a given simulation JSON.
 
     """
     translated_dict = None
@@ -537,7 +537,7 @@ def _translate_simulation_json(
         translated_dict = translation_func(input_params, mesh_unit)
     except Flow360TranslationError as err:
         raise ValueError(str(err)) from err
-    except Exception as err:  # tranlsation itself is not supposed to raise any other exception
+    except Exception as err:  # translation itself is not supposed to raise any other exception
         raise ValueError(
             f"Unexpected error translating to {target_name} json: " + str(err)
         ) from err
diff --git a/flow360/component/simulation/translator/solver_translator.py b/flow360/component/simulation/translator/solver_translator.py
index 137c9060f..783ba68d9 100644
--- a/flow360/component/simulation/translator/solver_translator.py
+++ b/flow360/component/simulation/translator/solver_translator.py
@@ -81,7 +81,18 @@
     translate_setting_and_apply_to_all_entities,
     update_dict_recursively,
 )
-from flow360.component.simulation.unit_system import LengthType
+from flow360.component.simulation.unit_system import (
+    CGS_unit_system,
+    LengthType,
+    SI_unit_system,
+    UnitSystem,
+    imperial_unit_system,
+    unit_system_manager,
+)
+from flow360.component.simulation.user_code import (
+    ValueOrExpression,
+    _solver_variable_values,
+)
 from flow360.component.simulation.utils import (
     is_exact_instance,
     is_instance_of_type_in_union,
@@ -219,6 +230,42 @@ def rotation_translator(model: Rotation):
     return volume_zone
 
 
+def _expression_has_single_solver_variable(expression: ValueOrExpression):
+    """
+    Check if the expression has a single solver variable.
+
+    This is used to determine if the expression is meant to dump a single native variable with potential requested target unit system. Instead of being a UDF.
+    """
+    return len(expression.solver_variable_names()) == 1
+
+
+def _get_expression_udf_name(output_field: ValueOrExpression):
+    if output_field.output_unit_system is None:
+        output_field.output_unit_system = unit_system_manager.current.name
+
+    # Note: 2 Possibilities:
+    # 1. `ValueOrExpression` Single Solver Variable:    Normal output with potential requested target unit system.
+    # 2. `ValueOrExpression` with UserVariable:         UDF written by user.
+    if _expression_has_single_solver_variable(output_field):
+        # 1.
+        field_name = output_field.solver_variable_names()[0]
+    else:
+        # 2.
+        raise NotImplementedError("UDF not supported yet")
+    # TODO: Better heuristic for output field name?
+    return field_name.partition("solution.")[2] + f"_{output_field.output_unit_system}"
+
+
+def _get_output_field_name(output_field: Union[ValueOrExpression, str]):
+    """Convert output field to string"""
+    if isinstance(output_field, str):
+        return output_field
+    if isinstance(output_field, ValueOrExpression):
+        return _get_expression_udf_name(output_field)
+
+    raise NotImplementedError()
+
+
 def translate_output_fields(
     output_model: Union[
         SurfaceOutput,
@@ -233,11 +280,15 @@ def translate_output_fields(
     ],
 ):
     """Get output fields"""
-    return {"outputFields": output_model.output_fields.items}
+    output_fields = []
+    for output_field in output_model.output_fields.items:
+        output_fields.append(_get_output_field_name(output_field))
+
+    return {"outputFields": output_fields}
 
 
 def surface_probe_setting_translation_func(entity: SurfaceProbeOutput):
-    """Translate non-entitties part of SurfaceProbeOutput"""
+    """Translate non-entities part of SurfaceProbeOutput"""
     dict_with_merged_output_fields = monitor_translator(entity)
     dict_with_merged_output_fields["surfacePatches"] = [
         surface.full_name for surface in entity.target_surfaces.stored_entities
@@ -351,11 +402,16 @@ def translate_volume_output(
         is_average=volume_output_class is TimeAverageVolumeOutput,
     )
     # Get outputFields
+    output_fields = []
+
+    for output_field in get_global_setting_from_first_instance(
+        output_params, volume_output_class, "output_fields"
+    ).items:
+        output_fields.append(_get_output_field_name(output_field))
+
     volume_output.update(
         {
-            "outputFields": get_global_setting_from_first_instance(
-                output_params, volume_output_class, "output_fields"
-            ).model_dump()["items"],
+            "outputFields": output_fields,
         }
     )
     return volume_output
@@ -506,6 +562,73 @@ def translate_acoustic_output(output_params: list):
     return None
 
 
+def get_solver_variable_with_converted_unit(
+    expression: ValueOrExpression, input_params: SimulationParams
+) -> ValueOrExpression:
+
+    def _get_dimensionalization_factor(
+        dimensionality, unit_system: UnitSystem, params: SimulationParams
+    ):
+        """
+        Get the dimensionalization factor for each component.
+        """
+        return params.convert_unit(
+            unit_system[dimensionality], target_system="flow360_v2"
+        ).value.item()
+
+    unit_system_name = (
+        expression.output_unit_system
+        if expression.output_unit_system
+        else unit_system_manager.current.name
+    )
+
+    solver_variable_name = expression.solver_variable_names()[0]
+    dimensionality = None
+
+    assert _expression_has_single_solver_variable(expression)
+
+    # Step 0: Figure out the dimensionality of the solver variable.
+    for variable_name, value in _solver_variable_values.items():
+        if variable_name == solver_variable_name:
+            dimensionality = value.units.dimensions
+            break
+
+    if not dimensionality:
+        raise Flow360TranslationError("Solver variable not found:", solver_variable_name)
+
+    # Step 1: Get the target unit for the solver variable.
+    unit_system = None
+
+    if unit_system_name == "SI":
+        unit_system = SI_unit_system
+    elif unit_system_name == "Imperial":
+        unit_system = imperial_unit_system
+    elif unit_system_name == "CGS":
+        unit_system = CGS_unit_system
+    else:
+        raise NotImplementedError()
+
+    # Step 2: Get the dimensionalization factor for each component.
+    conversion_constant = _get_dimensionalization_factor(dimensionality, unit_system, input_params)
+    return expression / conversion_constant
+
+
+def assemble_udf_expression(expression: ValueOrExpression, input_params: SimulationParams):
+    solver_variable_name = expression.solver_variable_names()[0]
+    for variable_name, value in _solver_variable_values.items():
+        if variable_name == solver_variable_name:
+            n_dim = value.value.ndim
+            break
+    if n_dim == 0:
+        return (
+            f"{_get_expression_udf_name(expression)} = {expression.to_solver_code(input_params)};"
+        )
+    assembled_expression = ""
+    for i_comp in range(len(value.value)):
+        assembled_expression += f"{_get_expression_udf_name(expression)}[{i_comp}] = ({expression.to_solver_code(input_params)})[{i_comp}];"
+    return assembled_expression
+
+
 def process_output_fields_for_udf(input_params: SimulationParams):
     """
     Process all output fields from different output types and generate additional
@@ -515,28 +638,39 @@ def process_output_fields_for_udf(input_params: SimulationParams):
         input_params: SimulationParams object containing outputs configuration
 
     Returns:
-        tuple: (all_field_names, generated_udfs) where:
-            - all_field_names is a set of all output field names
+        tuple: (all_fields, generated_udfs) where:
+            - all_fields is a set of all output field names
             - generated_udfs is a list of UserDefinedField objects for dimensioned fields
     """
 
     # Collect all output field names from all output types
-    all_field_names = set()
+    all_fields = set()
 
     if input_params.outputs:
         for output in input_params.outputs:
             if hasattr(output, "output_fields") and output.output_fields:
-                all_field_names.update(output.output_fields.items)
+                all_fields.update(output.output_fields.items)
 
     if isinstance(input_params.operating_condition, LiquidOperatingCondition):
-        all_field_names.add("velocity_magnitude")
+        all_fields.add("velocity_magnitude")
 
     # Generate UDFs for dimensioned fields
     generated_udfs = []
-    for field_name in all_field_names:
-        udf_expression = generate_predefined_udf(field_name, input_params)
-        if udf_expression:
-            generated_udfs.append(UserDefinedField(name=field_name, expression=udf_expression))
+    for field in all_fields:
+        if isinstance(field, str):
+            udf_expression = generate_predefined_udf(field, input_params)
+            if udf_expression:
+                generated_udfs.append(UserDefinedField(name=field, expression=udf_expression))
+
+        elif isinstance(field, ValueOrExpression):
+            converted_expression = get_solver_variable_with_converted_unit(field, input_params)
+            assembled_string = assemble_udf_expression(converted_expression, input_params)
+            generated_udfs.append(
+                UserDefinedField(
+                    name=_get_expression_udf_name(field),
+                    expression=assembled_string,
+                )
+            )
 
     return generated_udfs
 
@@ -979,7 +1113,7 @@ def get_solver_json(
     ##:: Step 1: Get geometry:
     if input_params.reference_geometry:
         geometry = inline_expressions_in_dict(
-            dump_dict(input_params.reference_geometry), input_params
+            remove_units_in_dict(dump_dict(input_params.reference_geometry)), input_params
         )
         geometry = remove_units_in_dict(geometry)
         translated["geometry"] = {}
@@ -1290,13 +1424,14 @@ def get_solver_json(
     )
 
     ##:: Step 4: Get outputs (has to be run after the boundaries are translated)
-
-    translated = translate_output(input_params, translated)
+    with input_params.unit_system:
+        translated = translate_output(input_params, translated)
 
     ##:: Step 5: Get user defined fields and auto-generated fields for dimensioned output
     translated["userDefinedFields"] = []
     # Add auto-generated UDFs for dimensioned fields
-    generated_udfs = process_output_fields_for_udf(input_params)
+    with input_params.unit_system:
+        generated_udfs = process_output_fields_for_udf(input_params)
 
     # Add user-specified UDFs and auto-generated UDFs for dimensioned fields
     for udf in [*input_params.user_defined_fields, *generated_udfs]:
diff --git a/flow360/component/simulation/translator/utils.py b/flow360/component/simulation/translator/utils.py
index 0f16496b8..b169f3bd2 100644
--- a/flow360/component/simulation/translator/utils.py
+++ b/flow360/component/simulation/translator/utils.py
@@ -119,10 +119,14 @@ def convert_tuples_to_lists(input_dict):
 
 def remove_units_in_dict(input_dict):
     """Remove units from a dimensioned value."""
+
+    def _is_unyt_or_unyt_like_obj(value):
+        return "value" in value.keys() and "units" in value.keys()
+
     unit_keys = {"value", "units"}
     if isinstance(input_dict, dict):
         new_dict = {}
-        if input_dict.keys() == unit_keys:
+        if _is_unyt_or_unyt_like_obj(input_dict):
             new_dict = input_dict["value"]
             if input_dict["units"].startswith("flow360_") is False:
                 raise ValueError(
@@ -130,7 +134,7 @@ def remove_units_in_dict(input_dict):
                 )
             return new_dict
         for key, value in input_dict.items():
-            if isinstance(value, dict) and value.keys() == unit_keys:
+            if isinstance(value, dict) and _is_unyt_or_unyt_like_obj(value):
                 if value["units"].startswith("flow360_") is False:
                     raise ValueError(
                         f"[Internal Error] Unit {value['units']} is not non-dimensionalized."
@@ -146,6 +150,7 @@ def remove_units_in_dict(input_dict):
 
 def inline_expressions_in_dict(input_dict, input_params):
     if isinstance(input_dict, dict):
+        print("Yes I am dict in the first place.")
         new_dict = {}
         if "expression" in input_dict.keys():
             expression = Expression(expression=input_dict["expression"])
diff --git a/flow360/component/simulation/unit_system.py b/flow360/component/simulation/unit_system.py
index f7c223504..0c6534a54 100644
--- a/flow360/component/simulation/unit_system.py
+++ b/flow360/component/simulation/unit_system.py
@@ -1596,7 +1596,17 @@ def defaults(self):
 
     def __getitem__(self, item):
         """to support [] access"""
-        return getattr(self, item)
+
+        try:
+            return getattr(self, item)
+        except TypeError:
+            # Allowing usage like [(mass)/(time)]
+            for attr_name, attr in vars(self).items():
+                if not isinstance(attr, unyt_quantity):
+                    continue
+                if attr.units.dimensions == item:
+                    return getattr(self, attr_name)
+        raise AttributeError(f"'{item}' is not a valid attribute of {self.__class__.__name__}. ")
 
     def system_repr(self):
         """(mass, length, time, temperature) string representation of the system"""
diff --git a/flow360/component/simulation/user_code.py b/flow360/component/simulation/user_code.py
index 2d0a4242a..2b395ec49 100644
--- a/flow360/component/simulation/user_code.py
+++ b/flow360/component/simulation/user_code.py
@@ -18,6 +18,7 @@
 _global_ctx: EvaluationContext = EvaluationContext(resolver)
 _user_variables: set[str] = set()
 _solver_variables: dict[str, str] = dict()
+_solver_variable_values: dict[str, unyt_quantity] = dict()
 
 
 def _is_number_string(s: str) -> bool:
@@ -83,12 +84,15 @@ def _convert_argument(value):
 
 
 class SerializedValueOrExpression(Flow360BaseModel):
-    type_name: Union[Literal["number"], Literal["expression"]] = pd.Field(None)
+    type_name: Literal["number", "expression"] = pd.Field()
     value: Optional[Union[Number, Iterable[Number]]] = pd.Field(None)
     units: Optional[str] = pd.Field(None)
     expression: Optional[str] = pd.Field(None)
     evaluated_value: Optional[Union[Number, Iterable[Number]]] = pd.Field(None)
     evaluated_units: Optional[str] = pd.Field(None)
+    output_unit_system: Optional[str] = pd.Field(
+        None, description="Results will be convert to this unit system if not None"
+    )
 
 
 # This is a wrapper to allow using ndarrays with pydantic models
@@ -254,11 +258,10 @@ def arctan(self):
 
 class UserVariable(Variable):
     @pd.model_validator(mode="after")
-    @classmethod
-    def update_context(cls, value):
-        _global_ctx.set(value.name, value.value)
-        _user_variables.add(value.name)
-        return value
+    def update_context(self):
+        _global_ctx.set(self.name, self.value)
+        _user_variables.add(self.name)
+        return self
 
     @pd.model_validator(mode="after")
     @classmethod
@@ -287,13 +290,20 @@ class SolverVariable(Variable):
     solver_name: Optional[str] = pd.Field(None)
 
     @pd.model_validator(mode="after")
-    @classmethod
-    def update_context(cls, value):
-        _global_ctx.set(value.name, value.value)
-        _solver_variables[value.name] = (
-            value.solver_name if value.solver_name is not None else value.name
+    def update_context(self):
+        _global_ctx.set(self.name, self.value)
+        _solver_variables[self.name] = (
+            self.solver_name if self.solver_name is not None else self.name
         )
-        return value
+        _solver_variable_values[self.name] = self.value
+        return self
+
+    def __hash__(self):
+        """
+        Support for set and deduplicate.
+        Can be removed if not used directly in output_fields.
+        """
+        return hash(self.model_dump_json())
 
 
 def _handle_syntax_error(se: SyntaxError, source: str):
@@ -321,20 +331,31 @@ def _handle_syntax_error(se: SyntaxError, source: str):
 
 class Expression(Flow360BaseModel, Evaluable):
     expression: str = pd.Field("")
-
-    model_config = pd.ConfigDict(validate_assignment=True)
+    output_unit_system: Optional[str] = pd.Field(None)
+    name: Optional[str] = pd.Field(
+        None, description="Name of the output variable if this is a output expression."
+    )
 
     @pd.model_validator(mode="before")
     @classmethod
     def _validate_expression(cls, value) -> Self:
+        name = None
+        output_unit_system = None
         if isinstance(value, str):
             expression = value
         elif isinstance(value, dict) and "expression" in value.keys():
             expression = value["expression"]
+            if "name" in value.keys():
+                name = value["name"]
+            if "output_unit_system" in value.keys():
+                output_unit_system = value["output_unit_system"]
         elif isinstance(value, Expression):
             expression = str(value)
+            name = value.name
+            output_unit_system = value.output_unit_system
         elif isinstance(value, Variable):
             expression = str(value)
+            name = value.name
         elif isinstance(value, np.ndarray) and not isinstance(value, unyt_array):
             if value.ndim == 0:
                 expression = str(value)
@@ -355,7 +376,7 @@ def _validate_expression(cls, value) -> Self:
             details = InitErrorDetails(type="value_error", ctx={"error": v_err})
             raise pd.ValidationError.from_exception_data("Expression value error", [details])
 
-        return {"expression": expression}
+        return {"expression": expression, "name": name, "output_unit_system": output_unit_system}
 
     def evaluate(
         self, context: EvaluationContext = None, strict: bool = True
@@ -380,6 +401,12 @@ def user_variable_names(self):
 
         return names
 
+    def solver_variable_names(self):
+        expr = expr_to_model(self.expression, _global_ctx)
+        names = expr.used_names()
+
+        return [name for name in names if name in _solver_variables]
+
     def to_solver_code(self, params):
         def translate_symbol(name):
             if name in _solver_variables:
@@ -563,7 +590,10 @@ def _deserialize(value) -> Self:
 
         def _serializer(value, info) -> dict:
             if isinstance(value, Expression):
-                serialized = SerializedValueOrExpression(type_name="expression")
+                serialized = SerializedValueOrExpression(
+                    type_name="expression",
+                    output_unit_system=value.output_unit_system,
+                )
 
                 serialized.expression = value.expression
 
diff --git a/flow360/component/simulation/utils.py b/flow360/component/simulation/utils.py
index e2bf9ceb2..9f687c347 100644
--- a/flow360/component/simulation/utils.py
+++ b/flow360/component/simulation/utils.py
@@ -57,34 +57,3 @@ def is_instance_of_type_in_union(obj, typ) -> bool:
 
     # Otherwise, do a normal isinstance check.
     return isinstance(obj, typ)
-
-
-class UnknownFloat(float):
-    def __new__(cls):
-        return super().__new__(cls, float("nan"))
-
-    def __repr__(self):
-        return "UnknownFloat()"
-
-    def __str__(self):
-        return "unknown"
-
-    def _return_unknown(self, *args):
-        return UnknownFloat()
-
-    __add__ = __radd__ = __sub__ = __rsub__ = _return_unknown
-    __mul__ = __rmul__ = __truediv__ = __rtruediv__ = _return_unknown
-    __floordiv__ = __rfloordiv__ = __mod__ = __rmod__ = _return_unknown
-    __pow__ = __rpow__ = _return_unknown
-    __neg__ = __pos__ = __abs__ = _return_unknown
-
-    def __eq__(self, other):
-        return False
-
-    def __ne__(self, other):
-        return True
-
-    __lt__ = __le__ = __gt__ = __ge__ = _return_unknown
-
-    def __bool__(self):
-        return False
diff --git a/flow360/component/simulation/validation/validation_output.py b/flow360/component/simulation/validation/validation_output.py
index 248606973..9b871c7a4 100644
--- a/flow360/component/simulation/validation/validation_output.py
+++ b/flow360/component/simulation/validation/validation_output.py
@@ -52,7 +52,7 @@ def extract_literal_values(annotation):
         allowed_items = natively_supported + additional_fields
 
         for item in output.output_fields.items:
-            if item not in allowed_items:
+            if item not in allowed_items and isinstance(item, str):
                 raise ValueError(
                     f"In `outputs`[{output_index}] {output.output_type}:, {item} is not a"
                     f" valid output field name. Allowed fields are {allowed_items}."
@@ -96,7 +96,7 @@ def _check_output_fields_valid_given_turbulence_model(params):
         if output.output_type in ("AeroAcousticOutput", "StreamlineOutput"):
             continue
         for item in output.output_fields.items:
-            if item in invalid_output_fields[turbulence_model]:
+            if isinstance(item, str) and item in invalid_output_fields[turbulence_model]:
                 raise ValueError(
                     f"In `outputs`[{output_index}] {output.output_type}:, {item} is not a valid"
                     f" output field when using turbulence model: {turbulence_model}."
diff --git a/flow360/component/simulation/variables/solution_variables.py b/flow360/component/simulation/variables/solution_variables.py
index 07237b7c8..8c72ec8e0 100644
--- a/flow360/component/simulation/variables/solution_variables.py
+++ b/flow360/component/simulation/variables/solution_variables.py
@@ -1,10 +1,18 @@
+from flow360.component.simulation import units as u
 from flow360.component.simulation.user_code import SolverVariable
 
-mut = SolverVariable(name="solution.mut", value=float("NaN"))  # Turbulent viscosity
-mu = SolverVariable(name="solution.mu", value=float("NaN"))  # Laminar viscosity
+mut = SolverVariable(
+    name="solution.mut", value=float("NaN") * u.kg / u.m / u.s, solver_name="mut"
+)  # Turbulent viscosity
+mu = SolverVariable(name="solution.mu", value=float("NaN") * u.kg / u.m / u.s)  # Laminar viscosity
+# TODO: So NaN will have to be serialized in the end now that we list them under the output_fields. @andrzej-krupka
+# TODO: We can bypass this problem by preprocessing `SolverVariables`
+# TODO: into `Expressions` and always serialize `Expression`
+# TODO: but I do not know if this is the right way to go since to_file_and_from_file test will fail.
 solutionNavierStokes = SolverVariable(
     name="solution.solutionNavierStokes", value=float("NaN")
 )  # Solution for N-S equation in conservative form
+# TODO: Do we need to support these? They have multiple dimensions inside a single SolverVariable
 residualNavierStokes = SolverVariable(
     name="solution.residualNavierStokes", value=float("NaN")
 )  # Residual for N-S equation in conservative form
diff --git a/tests/simulation/outputs/test_output_fields.py b/tests/simulation/outputs/test_output_fields.py
index 5ffbb7fbb..1f693e2ef 100644
--- a/tests/simulation/outputs/test_output_fields.py
+++ b/tests/simulation/outputs/test_output_fields.py
@@ -2,6 +2,7 @@
 
 import flow360 as fl
 from flow360 import SI_unit_system, u
+from flow360.component.simulation.framework.param_utils import AssetCache
 from flow360.component.simulation.outputs.output_fields import generate_predefined_udf
 
 
diff --git a/tests/simulation/test_expressions.py b/tests/simulation/test_expressions.py
index b8beaa4b4..84cba10ad 100644
--- a/tests/simulation/test_expressions.py
+++ b/tests/simulation/test_expressions.py
@@ -18,8 +18,9 @@
 from flow360.component.simulation.framework.base_model import Flow360BaseModel
 from flow360.component.simulation.framework.param_utils import AssetCache
 from flow360.component.simulation.models.material import Water, aluminum
-from flow360.component.simulation.outputs.outputs import SurfaceOutput
+from flow360.component.simulation.outputs.outputs import SurfaceOutput, VolumeOutput
 from flow360.component.simulation.primitives import GenericVolume, Surface
+from flow360.component.simulation.simulation_params import SimulationParams
 from flow360.component.simulation.unit_system import (
     AbsoluteTemperatureType,
     AngleType,
@@ -51,6 +52,7 @@
     UserVariable,
     ValueOrExpression,
 )
+from tests.utils import to_file_from_file_test
 
 
 @pytest.fixture(autouse=True)
@@ -667,6 +669,15 @@ def test_cyclic_dependencies():
         x.value = x
 
 
+def test_to_file_from_file_expression():
+    with SI_unit_system:
+        params = SimulationParams(
+            outputs=[VolumeOutput(output_fields=[solution.mut])],
+        )
+
+    to_file_from_file_test(params)
+
+
 def test_auto_alias():
     class TestModel(Flow360BaseModel):
         field: ValueOrExpression[VelocityType] = pd.Field()
diff --git a/tests/simulation/translator/ref/Flow360_volume_output_mut_dimensioned.json b/tests/simulation/translator/ref/Flow360_volume_output_mut_dimensioned.json
new file mode 100644
index 000000000..e640f9b2c
--- /dev/null
+++ b/tests/simulation/translator/ref/Flow360_volume_output_mut_dimensioned.json
@@ -0,0 +1,103 @@
+ {
+    "boundaries": {},
+    "freestream": {
+        "Mach": 0.029386353651012956,
+        "Temperature": 288.15,
+        "alphaAngle": 0.0,
+        "betaAngle": 0.0,
+        "muRef": 4.292321046986499e-08
+    },
+    "initialCondition": {
+        "p": "p",
+        "rho": "rho",
+        "type": "initialCondition",
+        "u": "u",
+        "v": "v",
+        "w": "w"
+    },
+    "navierStokesSolver": {
+        "CFLMultiplier": 1.0,
+        "absoluteTolerance": 1e-10,
+        "equationEvalFrequency": 1,
+        "kappaMUSCL": -1.0,
+        "limitPressureDensity": false,
+        "limitVelocity": false,
+        "linearSolver": {
+            "maxIterations": 30
+        },
+        "lowMachPreconditioner": false,
+        "maxForceJacUpdatePhysicalSteps": 0,
+        "modelType": "Compressible",
+        "numericalDissipationFactor": 1.0,
+        "orderOfAccuracy": 2,
+        "relativeTolerance": 0.0,
+        "updateJacobianFrequency": 4
+    },
+    "outputRescale": {
+        "velocityScale": 1.0
+    },
+    "timeStepping": {
+        "CFL": {
+            "convergenceLimitingFactor": 0.25,
+            "max": 10000.0,
+            "maxRelativeChange": 1.0,
+            "min": 0.1,
+            "type": "adaptive"
+        },
+        "maxPseudoSteps": 2000,
+        "orderOfAccuracy": 2,
+        "physicalSteps": 1,
+        "timeStepSize": "inf"
+    },
+    "turbulenceModelSolver": {
+        "CFLMultiplier": 2.0,
+        "DDES": false,
+        "ZDES": false,
+        "absoluteTolerance": 1e-08,
+        "equationEvalFrequency": 4,
+        "gridSizeForLES": "maxEdgeLength",
+        "linearSolver": {
+            "maxIterations": 20
+        },
+        "maxForceJacUpdatePhysicalSteps": 0,
+        "modelConstants": {
+            "C_DES": 0.72,
+            "C_cb1": 0.1355,
+            "C_cb2": 0.622,
+            "C_d": 8.0,
+            "C_min_rd": 10.0,
+            "C_sigma": 0.6666666666666666,
+            "C_t3": 1.2,
+            "C_t4": 0.5,
+            "C_v1": 7.1,
+            "C_vonKarman": 0.41,
+            "C_w2": 0.3
+        },
+        "modelType": "SpalartAllmaras",
+        "orderOfAccuracy": 2,
+        "quadraticConstitutiveRelation": false,
+        "reconstructionGradientLimiter": 0.5,
+        "relativeTolerance": 0.0,
+        "rotationCorrection": false,
+        "updateJacobianFrequency": 4
+    },
+    "userDefinedFields": [
+        {
+            "expression": "mut_SI = (whatever_mut_full_name * 0.0023988860123275884);",
+            "name": "mut_SI"
+        }
+    ],
+    "usingLiquidAsMaterial": false,
+    "volumeOutput": {
+        "animationFrequency": -1,
+        "animationFrequencyOffset": 0,
+        "animationFrequencyTimeAverage": -1,
+        "animationFrequencyTimeAverageOffset": 0,
+        "computeTimeAverages": false,
+        "outputFields": [
+            "mut_SI"
+        ],
+        "outputFormat": "paraview",
+        "startAverageIntegrationStep": -1
+    }
+}
\ No newline at end of file
diff --git a/tests/simulation/translator/test_solver_translator.py b/tests/simulation/translator/test_solver_translator.py
index eac407525..eb94a3597 100644
--- a/tests/simulation/translator/test_solver_translator.py
+++ b/tests/simulation/translator/test_solver_translator.py
@@ -108,10 +108,12 @@
 
 assertions = unittest.TestCase("__init__")
 
+from flow360.component.simulation.framework.param_utils import AssetCache
 from flow360.component.simulation.framework.updater_utils import compare_values
 from flow360.component.simulation.models.volume_models import AngleExpression, Rotation
 from flow360.component.simulation.primitives import GenericVolume
 from flow360.component.simulation.time_stepping.time_stepping import Unsteady
+from flow360.component.simulation.variables import solution_variables as solution
 
 
 @pytest.fixture()
@@ -612,3 +614,19 @@ def test_liquid_simulation_translation():
         # Flow360 time to seconds = 1m/(50m/s) = 0.02 s
         # t_seconds = (0.02 s * t)
     translate_and_compare(param, mesh_unit=1 * u.m, ref_json_file="Flow360_liquid_rotation_dd.json")
+
+
+def test_solution_variable_output():
+    with SI_unit_system:
+        params = SimulationParams(
+            operating_condition=AerospaceCondition(velocity_magnitude=10),
+            outputs=[VolumeOutput(output_fields=[solution.mut])],
+            private_attribute_asset_cache=AssetCache(project_length_unit="m"),
+        )
+
+    translate_and_compare(
+        params,
+        mesh_unit=1 * u.m,
+        debug=True,
+        ref_json_file="Flow360_volume_output_mut_dimensioned.json",
+    )