Add eulerian multiphase

CfdOF/CfdTools.py

@@ -288,6 +288,13 @@ def getScalarTransportFunctionsGroup(analysis_object):
     return group
 
 
+# Reactions
+def getReactionModels(analysis_object):
+    group = []
+    # TODO complete
+    return group
+
+
 # Mesh
 def getMeshRefinementObjs(mesh_obj):
     from CfdOF.Mesh.CfdMeshRefinement import CfdMeshRefinement

CfdOF/Solve/CfdCaseWriterFoam.py

@@ -45,6 +45,7 @@ def __init__(self, analysis_obj):
         self.bc_group = CfdTools.getCfdBoundaryGroup(analysis_obj)
         self.initial_conditions = CfdTools.getInitialConditions(analysis_obj)
         self.reporting_functions = CfdTools.getReportingFunctionsGroup(analysis_obj)
+        self.reaction_models = CfdTools.getReactionModels(analysis_obj)
         self.scalar_transport_objs = CfdTools.getScalarTransportFunctionsGroup(analysis_obj)
         self.porous_zone_objs = CfdTools.getPorousZoneObjects(analysis_obj)
         self.initialisation_zone_objs = CfdTools.getInitialisationZoneObjects(analysis_obj)
@@ -90,6 +91,7 @@ def writeCase(self):
             'boundaries': dict((b.Label, CfdTools.propsToDict(b)) for b in self.bc_group),
             'reportingFunctions': dict((fo.Label, CfdTools.propsToDict(fo)) for fo in self.reporting_functions),
             'reportingFunctionsEnabled': False,
+            'reactionsEnabled': False,
             'scalarTransportFunctions': dict((st.Label, CfdTools.propsToDict(st)) for st in self.scalar_transport_objs),
             'scalarTransportFunctionsEnabled': False,
             'dynamicMesh': {},
@@ -138,6 +140,10 @@ def writeCase(self):
             cfdMessage('Dynamic mesh adaptation rule present\n')
             self.processDynamicMeshRefinement()
 
+        if self.reaction_models:
+            cfdMessage('Reactions present\n')
+            self.processReactions()
+
         self.settings['createPatchesFromSnappyBaffles'] = False
         self.settings['createPatchesForPeriodics'] = False
         cfdMessage("Matching boundary conditions ...\n")
@@ -201,6 +207,11 @@ def getSolverName(self):
                     raise RuntimeError("Only isothermal analysis is currently supported for free surface flow simulation.")
             else:
                 raise RuntimeError("Only transient analysis is supported for free surface flow simulation.")
+        elif self.physics_model.Phase == 'Eulerian':
+            if len(self.material_objs) >= 2:
+                solver = 'multiphaseEulerFoam'
+            else:
+                raise RuntimeError("At least two material objects must be present for Eulerian multiphase simulation.")
         else:
             raise RuntimeError(self.physics_model.Phase + " phase model currently not supported.")
 
@@ -331,15 +342,15 @@ def processBoundaryConditions(self):
                 beta = (1-wireDiam/spacing)**2
                 bc['PressureDropCoeff'] = CD*(1-beta)
 
-            if settings['solver']['SolverName'] in ['interFoam', 'multiphaseInterFoam']:
+            if settings['solver']['SolverName'] in ['interFoam', 'multiphaseInterFoam', 'multiphaseEulerFoam']:
                 # Make sure the first n-1 alpha values exist, and write the n-th one
                 # consistently for multiphaseInterFoam
                 sum_alpha = 0.0
                 alphas_new = {}
                 for i, m in enumerate(settings['fluidProperties']):
                     alpha_name = m['Name']
                     if i == len(settings['fluidProperties']) - 1:
-                        if settings['solver']['SolverName'] == 'multiphaseInterFoam':
+                        if settings['solver']['SolverName'] == 'multiphaseInterFoam' or settings['solver']['SolverName'] == 'multiphaseEulerFoam':
                             alphas_new[alpha_name] = 1.0 - sum_alpha
                     else:
                         alpha = Units.Quantity(bc.get('VolumeFractions', {}).get(alpha_name, '0')).Value
@@ -417,15 +428,15 @@ def processInitialConditions(self):
         if settings['solver']['SolverName'] in ['simpleFoam', 'porousSimpleFoam', 'pimpleFoam', 'SRFSimpleFoam']:
             mat_prop = settings['fluidProperties'][0]
             initial_values['KinematicPressure'] = initial_values['Pressure'] / mat_prop['Density']
-        if settings['solver']['SolverName'] in ['interFoam', 'multiphaseInterFoam']:
+        if settings['solver']['SolverName'] in ['interFoam', 'multiphaseInterFoam', 'multiphaseEulerFoam']:
             # Make sure the first n-1 alpha values exist, and write the n-th one
             # consistently for multiphaseInterFoam
             sum_alpha = 0.0
             alphas_new = {}
             for i, m in enumerate(settings['fluidProperties']):
                 alpha_name = m['Name']
                 if i == len(settings['fluidProperties'])-1:
-                    if settings['solver']['SolverName'] == 'multiphaseInterFoam':
+                    if settings['solver']['SolverName'] == 'multiphaseInterFoam' or settings['solver']['SolverName'] == 'multiphaseEulerFoam':
                         alphas_new[alpha_name] = 1.0-sum_alpha
                 else:
                     alpha = Units.Quantity(initial_values.get('VolumeFractions', {}).get(alpha_name, '0')).Value
@@ -570,7 +581,7 @@ def processReportingFunctions(self):
                 rf['Pitch'] = tuple(p for p in pitch_axis)
 
             settings['reportingFunctions'][name]['ProbePosition'] = tuple(
-                Units.Quantity(p, Units.Length).getValueAs('m') 
+                Units.Quantity(p, Units.Length).getValueAs('m')
                 for p in settings['reportingFunctions'][name]['ProbePosition'])
 
     def processScalarTransportFunctions(self):
@@ -584,6 +595,11 @@ def processScalarTransportFunctions(self):
             stf['InjectionPoint'] = tuple(
                 Units.Quantity(p, Units.Length).getValueAs('m') for p in stf['InjectionPoint'])
 
+    # Process reactions and reactions models
+    def processReactions(self):
+        settings = self.settings
+        settings['reactionsEnabled'] = True
+
     # Mesh related
     def processDynamicMeshRefinement(self):
         settings = self.settings
@@ -592,14 +608,14 @@ def processDynamicMeshRefinement(self):
         # Check whether cellLevel supported
         if self.mesh_obj.MeshUtility not in ['cfMesh', 'snappyHexMesh']:
             raise RuntimeError("Dynamic mesh refinement is only supported by cfMesh and snappyHexMesh")
-    
+
         # Check whether 2D extrusion present
         mesh_refinements = CfdTools.getMeshRefinementObjs(self.mesh_obj)
         for mr in mesh_refinements:
             if mr.Extrusion:
                 if mr.ExtrusionType == '2DPlanar' or mr.ExtrusionType == '2DWedge':
                     raise RuntimeError("Dynamic mesh refinement will not work with 2D or wedge mesh")
-        
+
         settings['dynamicMesh'] = CfdTools.propsToDict(self.dynamic_mesh_refinement_obj)
         if isinstance(self.dynamic_mesh_refinement_obj.Proxy, CfdDynamicMeshRefinement.CfdDynamicMeshShockRefinement):
             settings['dynamicMesh']['Type'] = 'shock'
@@ -669,7 +685,7 @@ def processPorousZoneProperties(self):
 
     def processInitialisationZoneProperties(self):
         settings = self.settings
-        
+
         for zone_name in settings['initialisationZones']:
             z = settings['initialisationZones'][zone_name]
 
@@ -684,7 +700,7 @@ def processInitialisationZoneProperties(self):
             if not z['TemperatureSpecified']:
                 del z['Temperature']
 
-            if settings['solver']['SolverName'] in ['interFoam', 'multiphaseInterFoam']:
+            if settings['solver']['SolverName'] in ['interFoam', 'multiphaseInterFoam', 'multiphaseEulerFoam']:
                 # Make sure the first n-1 alpha values exist, and write the n-th one
                 # consistently for multiphaseInterFoam
                 sum_alpha = 0.0
@@ -693,14 +709,14 @@ def processInitialisationZoneProperties(self):
                     for i, m in enumerate(settings['fluidProperties']):
                         alpha_name = m['Name']
                         if i == len(settings['fluidProperties'])-1:
-                            if settings['solver']['SolverName'] == 'multiphaseInterFoam':
+                            if settings['solver']['SolverName'] == 'multiphaseInterFoam' or settings['solver']['SolverName'] == 'multiphaseEulerFoam':
                                 alphas_new[alpha_name] = 1.0-sum_alpha
                         else:
                             alpha = Units.Quantity(z['VolumeFractions'].get(alpha_name, '0')).Value
                             alphas_new[alpha_name] = alpha
                             sum_alpha += alpha
                     z['VolumeFractions'] = alphas_new
-        
+
             if settings['solver']['SolverName'] in ['simpleFoam', 'porousSimpleFoam', 'pimpleFoam', 'SRFSimpleFoam']:
                 if 'Pressure' in z:
                     z['KinematicPressure'] = z['Pressure']/settings['fluidProperties'][0]['Density']

CfdOF/Solve/CfdPhasePhysicsSelection.py

@@ -0,0 +1,231 @@
+# ***************************************************************************
+# *                                                                         *
+# *   Copyright (c) 2024 Jonathan Bergh <[email protected]>          *
+# *                                                                         *
+# *   This program is free software: you can redistribute it and/or modify  *
+# *   it under the terms of the GNU Lesser General Public License as        *
+# *   published by the Free Software Foundation, either version 3 of the    *
+# *   License, or (at your option) any later version.                       *
+# *                                                                         *
+# *   This program is distributed in the hope that it will be useful,       *
+# *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+# *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+# *   GNU Lesser General Public License for more details.                   *
+# *                                                                         *
+# *   You should have received a copy of the GNU Lesser General Public      *
+# *   License along with this program.  If not,                             *
+# *   see <https://www.gnu.org/licenses/>.                                  *
+# *                                                                         *
+# ***************************************************************************
+
+
+import os
+import os.path
+import FreeCAD
+if FreeCAD.GuiUp:
+    import FreeCADGui
+    from PySide import QtCore
+from CfdOF import CfdTools
+from CfdOF.CfdTools import addObjectProperty
+
+
+def makeCfdPhasePhysicsSelection(name="PhasePhysicsModel"):
+    # DocumentObjectGroupPython, FeaturePython, GeometryPython
+    obj = FreeCAD.ActiveDocument.addObject("Part::FeaturePython", name)
+    CfdPhasePhysicsModel(obj)
+    if FreeCAD.GuiUp:
+        ViewProviderCfdPhasePhysicsSelection(obj.ViewObject)
+    return obj
+
+
+class CommandCfdPhasePhysicsSelection:
+    """ CFD physics selection command definition """
+
+    def GetResources(self):
+        icon_path = os.path.join(CfdTools.getModulePath(), "Gui", "Icons", "phasephysics.svg")
+        return {'Pixmap': icon_path,
+                'MenuText': QtCore.QT_TRANSLATE_NOOP("Cfd_PhasePhysicsModel", "Select Eulerian phase models"),
+                'Accel': "",
+                'ToolTip': QtCore.QT_TRANSLATE_NOOP("Cfd_PhysicsModel", "Select the Eulerian multiphase physics models")}
+
+    def IsActive(self):
+        return CfdTools.getActiveAnalysis() is not None
+
+    def Activated(self):
+        FreeCAD.ActiveDocument.openTransaction("Choose appropriate phase physics model")
+        is_present = False
+        members = CfdTools.getActiveAnalysis().Group
+        for i in members:
+            if isinstance(i.Proxy, CfdPhasePhysicsModel):
+                FreeCADGui.activeDocument().setEdit(i.Name)
+                is_present = True
+
+        # Allow to re-create if deleted
+        if not is_present:
+            FreeCADGui.doCommand("")
+            FreeCADGui.doCommand("from CfdOF.Solve import CfdPhasePhysicsSelection")
+            FreeCADGui.doCommand("from CfdOF import CfdTools")
+            FreeCADGui.doCommand(
+                "CfdTools.getActiveAnalysis().addObject(CfdPhasePhysicsSelection.makeCfdPhasePhysicsSelection())")
+            FreeCADGui.ActiveDocument.setEdit(FreeCAD.ActiveDocument.ActiveObject.Name)
+
+
+class CfdPhasePhysicsModel:
+    """ The Eulerian multiphase specific physics models """
+    def __init__(self, obj):
+        obj.Proxy = self
+        self.Type = "PhasePhysicsModel"
+        self.initProperties(obj)
+
+    def initProperties(self, obj):
+
+        if addObjectProperty(obj, "Drag", ['SchillerNeuman'], "App::PropertyEnumeration", "Phase Physics modelling",
+                             "Drag force model"):
+            obj.Drag = 'SchillerNeuman'
+
+        if addObjectProperty(obj, "Lift", [''],
+                             "App::PropertyEnumeration", "Phase Physics modelling", "Lift force model"):
+            obj.Lift = ''
+
+        if addObjectProperty(obj, "SurfaceTension", ['Constant'], "App::PropertyEnumeration", "Phase Physics modelling",
+                             "Surface tension force model"):
+            obj.SurfaceTension = 'Constant'
+
+        if addObjectProperty(obj, "VirtualMass", ['ConstantCoefficient'],
+                             "App::PropertyEnumeration", "Phase Physics modelling", "Virtual mass force model"):
+            obj.VirtualMass = 'ConstantCoefficient'
+
+        if addObjectProperty(obj, "HeatTransfer", ['RanzMarshall'],
+                             "App::PropertyEnumeration", "Phase Physics modelling", "Interphase heat transfer model"):
+            obj.HeatTransfer = 'RanzMarshall'
+
+        if addObjectProperty(obj, "PhaseTransfer", [''],
+                             "App::PropertyEnumeration", "Phase Physics modelling", "Interphase mass transfer model"):
+            obj.PhaseTransfer = ''
+
+        if addObjectProperty(obj, "WallLubrication", [''],
+                             "App::PropertyEnumeration", "Phase Physics modelling", "Wall lubrication model"):
+            obj.WallLubrication = ''
+
+        if addObjectProperty(obj, "TurbulentDispersion", [''],
+                             "App::PropertyEnumeration", "Phase Physics modelling", "Turbulent dispersion model"):
+            obj.TurbulentDispersion = ''
+
+        if addObjectProperty(obj, "InterfaceCompression", [''],
+                             "App::PropertyEnumeration", "Phase Physics modelling", "Interface compression model"):
+            obj.InterfaceCompression = ''
+
+        # SRF model
+        addObjectProperty(obj, 'SRFModelEnabled', False, "App::PropertyBool", "Reference frame",
+                          "Single Rotating Frame model enabled")
+
+        addObjectProperty(obj, 'SRFModelRPM', '0', "App::PropertyQuantity", "Reference frame", "Rotational speed")
+
+        addObjectProperty(obj, 'SRFModelCoR', FreeCAD.Vector(0, 0, 0), "App::PropertyPosition", "Reference frame",
+                          "Centre of rotation (SRF)")
+
+        addObjectProperty(obj, 'SRFModelAxis', FreeCAD.Vector(0, 0, 0), "App::PropertyPosition", "Reference frame",
+                          "Axis of rotation (SRF)")
+
+    def onDocumentRestored(self, obj):
+        self.initProperties(obj)
+
+
+class _CfdPhasePhysicsModel:
+    """ Backward compatibility for old class name when loading from file """
+    def onDocumentRestored(self, obj):
+        CfdPhasePhysicsModel(obj)
+
+    def __getstate__(self):
+        return None
+
+    def __setstate__(self, state):
+        return None
+
+    # dumps and loads replace __getstate__ and __setstate__ post v. 0.21.2
+    def dumps(self):
+        return None
+
+    def loads(self, state):
+        return None
+
+
+class ViewProviderCfdPhasePhysicsSelection:
+    def __init__(self, vobj):
+        vobj.Proxy = self
+        self.taskd = None
+
+    def getIcon(self):
+        icon_path = os.path.join(CfdTools.getModulePath(), "Gui", "Icons", "phasephysics.svg")
+        return icon_path
+
+    def attach(self, vobj):
+        self.ViewObject = vobj
+        self.Object = vobj.Object
+        self.bubbles = None
+
+    def updateData(self, obj, prop):
+        analysis_obj = CfdTools.getParentAnalysisObject(obj)
+        if analysis_obj and not analysis_obj.Proxy.loading:
+            analysis_obj.NeedsCaseRewrite = True
+
+    def onChanged(self, vobj, prop):
+        return
+
+    def setEdit(self, vobj, mode):
+        from CfdOF.Solve import TaskPanelCfdPhasePhysicsSelection
+        import importlib
+        importlib.reload(TaskPanelCfdPhasePhysicsSelection)
+        self.taskd = TaskPanelCfdPhasePhysicsSelection.TaskPanelCfdPhasePhysicsSelection(self.Object)
+        self.taskd.obj = vobj.Object
+        FreeCADGui.Control.showDialog(self.taskd)
+        return True
+
+    def doubleClicked(self, vobj):
+        # Make sure no other task dialog still active
+        doc = FreeCADGui.getDocument(vobj.Object.Document)
+        if not doc.getInEdit():
+            doc.setEdit(vobj.Object.Name)
+        else:
+            FreeCAD.Console.PrintError('Task dialog already active\n')
+            FreeCADGui.Control.showTaskView()
+        return True
+
+    def unsetEdit(self, vobj, mode):
+        if self.taskd:
+            self.taskd.closing()
+            self.taskd = None
+        FreeCADGui.Control.closeDialog()
+
+    def __getstate__(self):
+        return None
+
+    def __setstate__(self, state):
+        return None
+
+    # dumps and loads replace __getstate__ and __setstate__ post v. 0.21.2
+    def dumps(self):
+        return None
+
+    def loads(self, state):
+        return None
+
+
+class _ViewProviderPhasePhysicsSelection:
+    """ Backward compatibility for old class name when loading from file """
+    def attach(self, vobj):
+        new_proxy = ViewProviderCfdPhasePhysicsSelection(vobj)
+        new_proxy.attach(vobj)
+
+    def __getstate__(self):
+        return None
+
+    def __setstate__(self, state):
+        return None
+
+    # dumps and loads replace __getstate__ and __setstate__ post v. 0.21.2
+    def dumps(self):
+        return None
+
+    def loads(self, state):
+        return None