From cf293180633b1a520a2568cf32f6eedd043cb050 Mon Sep 17 00:00:00 2001 From: Uwe Date: Fri, 17 Mar 2023 17:30:21 +0100 Subject: [PATCH] [FEM] add an example for a turbulent flow --- src/Mod/Fem/CMakeLists.txt | 1 + .../equation_flow_turbulent_elmer_2D.py | 293 ++++++++++++++++++ src/Mod/Fem/femexamples/manager.py | 2 + 3 files changed, 296 insertions(+) create mode 100644 src/Mod/Fem/femexamples/equation_flow_turbulent_elmer_2D.py diff --git a/src/Mod/Fem/CMakeLists.txt b/src/Mod/Fem/CMakeLists.txt index 233939176b..e327feb1cc 100755 --- a/src/Mod/Fem/CMakeLists.txt +++ b/src/Mod/Fem/CMakeLists.txt @@ -78,6 +78,7 @@ SET(FemExamples_SRCS femexamples/equation_electrostatics_capacitance_two_balls.py femexamples/equation_electrostatics_electricforce_elmer_nongui6.py femexamples/equation_flow_elmer_2D.py + femexamples/equation_flow_turbulent_elmer_2D.py femexamples/equation_flux_elmer.py femexamples/equation_magnetodynamics_elmer.py femexamples/equation_magnetodynamics_2D_elmer.py diff --git a/src/Mod/Fem/femexamples/equation_flow_turbulent_elmer_2D.py b/src/Mod/Fem/femexamples/equation_flow_turbulent_elmer_2D.py new file mode 100644 index 0000000000..e21510bc75 --- /dev/null +++ b/src/Mod/Fem/femexamples/equation_flow_turbulent_elmer_2D.py @@ -0,0 +1,293 @@ +# *************************************************************************** +# * Copyright (c) 2023 Uwe Stöhr * +# * * +# * This file is part of the FreeCAD CAx development system. * +# * * +# * This program is free software; you can redistribute it and/or modify * +# * it under the terms of the GNU Lesser General Public License (LGPL) * +# * as published by the Free Software Foundation; either version 2 of * +# * the License, or (at your option) any later version. * +# * for detail see the LICENCE text file. * +# * * +# * 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 Library General Public License for more details. * +# * * +# * You should have received a copy of the GNU Library General Public * +# * License along with this program; if not, write to the Free Software * +# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 * +# * USA * +# * * +# *************************************************************************** + +import sys +import FreeCAD +from FreeCAD import Placement +from FreeCAD import Rotation +from FreeCAD import Vector + +import Draft +import ObjectsFem + +from BOPTools import SplitFeatures +from . import manager +from .manager import get_meshname +from .manager import init_doc + +def get_information(): + return { + "name": "Turbulent Flow - Elmer 2D", + "meshtype": "solid", + "meshelement": "Tet10", + "constraints": ["initial pressure", "initial temperature", "initial velocity", + "temperature", "velocity"], + "solvers": ["elmer"], + "material": "fluid", + "equations": ["flow", "heat"] + } + +def get_explanation(header=""): + return header + """ + +To run the example from Python console use: +from femexamples.equation_flow_turbulent_elmer_2D import setup +setup() + +Flow and Heat equation - Elmer solver + +""" + +def setup(doc=None, solvertype="elmer"): + + # init FreeCAD document + if doc is None: + doc = init_doc() + + # explanation object + # just keep the following line and change text string in get_explanation method + manager.add_explanation_obj(doc, get_explanation(manager.get_header(get_information()))) + + # geometric objects + + # the wire defining the pipe volume in 2D + p1 = Vector(400, 0, -50.000) + p2 = Vector(400, 0, -150.000) + p3 = Vector(1200, 0, -150.000) + p4 = Vector(1200, 0, 50.000) + p5 = Vector(0, 0, 50.000) + p6 = Vector(0, 0, -50.000) + wire = Draft.make_wire([p1, p2, p3, p4, p5, p6], closed=True) + wire.Label = "Wire" + + # the circle defining the heating rod + pCirc = Vector(160, 0, 0) + axisCirc = Vector(1, 0, 0) + placementCircle = Placement(pCirc, Rotation(axisCirc, 90)) + circle = Draft.make_circle(10, placement=placementCircle) + circle.Label = "HeatingRod" + circle.ViewObject.Visibility = False + + # a link of the circle + circleLink = doc.addObject("App::Link", "Link-HeatingRod") + circleLink.LinkTransform = True + circleLink.LinkedObject = circle + + # cut rod from wire to get volume of fluid + cut = doc.addObject("Part::Cut", "Cut") + cut.Base = wire + cut.Tool = circleLink + cut.ViewObject.Visibility = False + + # BooleanFregments object to combine cut with rod + BooleanFragments = SplitFeatures.makeBooleanFragments(name="BooleanFragments") + BooleanFragments.Objects = [cut, circle] + + # set view + doc.recompute() + if FreeCAD.GuiUp: + BooleanFragments.ViewObject.Transparency = 50 + BooleanFragments.ViewObject.Document.activeView().viewFront() + BooleanFragments.ViewObject.Document.activeView().fitAll() + + # analysis + analysis = ObjectsFem.makeAnalysis(doc, "Analysis") + if FreeCAD.GuiUp: + import FemGui + FemGui.setActiveAnalysis(analysis) + + # solver + if solvertype == "elmer": + solver_obj = ObjectsFem.makeSolverElmer(doc, "SolverElmer") + equation_flow = ObjectsFem.makeEquationFlow(doc, solver_obj) + equation_heat = ObjectsFem.makeEquationHeat(doc, solver_obj) + else: + FreeCAD.Console.PrintWarning( + "Unknown or unsupported solver type: {}. " + "No solver object was created.\n".format(solvertype) + ) + return doc + analysis.addObject(solver_obj) + + # solver settings + equation_flow.IdrsParameter = 3 + equation_flow.LinearIterations = 250 + equation_flow.LinearIterativeMethod = "Idrs" + equation_flow.LinearPreconditioning = "ILU1" + equation_flow.setExpression("LinearTolerance", "1e-6") + equation_flow.NonlinearIterations = 30 + equation_flow.NonlinearNewtonAfterIterations = 30 + equation_flow.setExpression("NonlinearTolerance", "1e-4") + equation_flow.RelaxationFactor = 0.1 + equation_heat.Convection = "Computed" + equation_heat.IdrsParameter = 3 + equation_heat.LinearIterations = 250 + equation_heat.LinearIterativeMethod = "Idrs" + equation_heat.LinearPreconditioning = "ILU1" + equation_heat.setExpression("LinearTolerance", "1e-6") + equation_heat.NonlinearIterations = 30 + equation_heat.NonlinearNewtonAfterIterations = 30 + equation_heat.setExpression("NonlinearTolerance", "1e-4") + equation_heat.Priority = 5 + equation_heat.RelaxationFactor = 0.1 + equation_heat.Stabilize = True + + # material + + # fluid + material_obj = ObjectsFem.makeMaterialFluid(doc, "Material_Fluid") + mat = material_obj.Material + mat["Name"] = "Water" + mat["Density"] = "998 kg/m^3" + mat["DynamicViscosity"] = "1.003e-3 kg/m/s" + mat["ThermalConductivity"] = "0.591 W/m/K" + mat["ThermalExpansionCoefficient"] = "2.07e-4 m/m/K" + mat["SpecificHeat"] = "4182 J/kg/K" + material_obj.Material = mat + material_obj.References = [(BooleanFragments, "Face2")] + analysis.addObject(material_obj) + + # tube wall + material_obj = ObjectsFem.makeMaterialSolid(doc, "Material_Wall") + mat = material_obj.Material + mat["Name"] = "Aluminum Generic" + mat["Density"] = "2700 kg/m^3" + mat["PoissonRatio"] = "0.35" + mat["ShearModulus"] = "25.0 GPa" + mat["UltimateTensileStrength"] = "310 MPa" + mat["YoungsModulus"] = "70000 MPa" + mat["ThermalConductivity"] = "237.0 W/m/K" + mat["ThermalExpansionCoefficient"] = "23.1 µm/m/K" + mat["SpecificHeat"] = "897.0 J/kg/K" + material_obj.Material = mat + material_obj.References = [(BooleanFragments, "Face1")] + analysis.addObject(material_obj) + + # constraint inlet velocity + FlowVelocity_Inlet = ObjectsFem.makeConstraintFlowVelocity(doc, "FlowVelocity_Inlet") + FlowVelocity_Inlet.References = [(BooleanFragments, "Edge5")] + FlowVelocity_Inlet.NormalDirection = Vector(-1, 0, 0) + FlowVelocity_Inlet.VelocityX = 0.020 + FlowVelocity_Inlet.VelocityXEnabled = True + FlowVelocity_Inlet.VelocityYEnabled = True + FlowVelocity_Inlet.VelocityZEnabled = True + analysis.addObject(FlowVelocity_Inlet) + + # constraint outlet velocity + FlowVelocity_Outlet = ObjectsFem.makeConstraintFlowVelocity(doc, "FlowVelocity_Outlet") + FlowVelocity_Outlet.References = [(BooleanFragments, "Edge6")] + FlowVelocity_Outlet.NormalDirection = Vector(1, 0, 0) + FlowVelocity_Outlet.VelocityYEnabled = True + FlowVelocity_Outlet.VelocityZEnabled = True + analysis.addObject(FlowVelocity_Outlet) + + # constraint wall velocity + FlowVelocity_Wall = ObjectsFem.makeConstraintFlowVelocity(doc, "FlowVelocity_Wall") + FlowVelocity_Wall.References = [ + (BooleanFragments, "Edge2"), + (BooleanFragments, "Edge3"), + (BooleanFragments, "Edge4"), + (BooleanFragments, "Edge7")] + FlowVelocity_Wall.NormalDirection = Vector(0, 0, -1) + FlowVelocity_Wall.VelocityXEnabled = True + FlowVelocity_Wall.VelocityYEnabled = True + FlowVelocity_Wall.VelocityZEnabled = True + analysis.addObject(FlowVelocity_Wall) + + # constraint initial velocity + FlowVelocity_Initial = ObjectsFem.makeConstraintInitialFlowVelocity(doc, "FlowVelocity_Initial") + FlowVelocity_Initial.References = [(BooleanFragments, "Face2")] + FlowVelocity_Initial.NormalDirection = Vector(0, -1, 0) + FlowVelocity_Initial.VelocityXEnabled = True + FlowVelocity_Initial.VelocityYEnabled = True + FlowVelocity_Initial.VelocityZEnabled = True + analysis.addObject(FlowVelocity_Initial) + + # constraint initial temperature + Temperature_Initial = ObjectsFem.makeConstraintInitialTemperature(doc, "Temperature_Initial") + Temperature_Initial.initialTemperature = 300.0 + analysis.addObject(Temperature_Initial) + + # constraint wall temperature + Temperature_Wall = ObjectsFem.makeConstraintTemperature(doc, "Temperature_Wall") + Temperature_Wall.Temperature = 300.0 + Temperature_Wall.NormalDirection = Vector(0, 0, -1) + Temperature_Wall.References = [ + (BooleanFragments, "Edge2"), + (BooleanFragments, "Edge3"), + (BooleanFragments, "Edge4"), + (BooleanFragments, "Edge7")] + analysis.addObject(Temperature_Wall) + + # constraint inlet temperature + Temperature_Inlet = ObjectsFem.makeConstraintTemperature(doc, "Temperature_Inlet") + Temperature_Inlet.Temperature = 350.0 + Temperature_Inlet.NormalDirection = Vector(-1, 0, 0) + Temperature_Inlet.References = [(BooleanFragments, "Edge5")] + analysis.addObject(Temperature_Inlet) + + # constraint heating rod temperature + Temperature_HeatingRod = ObjectsFem.makeConstraintTemperature(doc, "Temperature_HeatingRod") + Temperature_HeatingRod.Temperature = 373.0 + Temperature_HeatingRod.NormalDirection = Vector(0, -1, 0) + Temperature_HeatingRod.References = [(BooleanFragments, "Edge1")] + analysis.addObject(Temperature_HeatingRod) + + # constraint initial pressure + Pressure_Initial = ObjectsFem.makeConstraintInitialPressure(doc, "Pressure_Initial") + Pressure_Initial.Pressure = "100.0 kPa" + Pressure_Initial.NormalDirection = Vector(0, -1, 0) + Pressure_Initial.References = [(BooleanFragments, "Face2")] + analysis.addObject(Pressure_Initial) + + # mesh + femmesh_obj = analysis.addObject(ObjectsFem.makeMeshGmsh(doc, get_meshname()))[0] + femmesh_obj.Part = BooleanFragments + femmesh_obj.ElementOrder = "1st" + femmesh_obj.CharacteristicLengthMax = "4 mm" + femmesh_obj.ViewObject.Visibility = False + + # mesh_region + mesh_region = ObjectsFem.makeMeshRegion(doc, femmesh_obj, name="MeshRegion") + mesh_region.CharacteristicLength = "2 mm" + mesh_region.References = [ + (BooleanFragments, "Edge1"), + (BooleanFragments, "Vertex2"), + (BooleanFragments, "Vertex4"), + (BooleanFragments, "Vertex6")] + mesh_region.ViewObject.Visibility = False + + # generate the mesh + from femmesh import gmshtools + gmsh_mesh = gmshtools.GmshTools(femmesh_obj, analysis) + try: + error = gmsh_mesh.create_mesh() + except Exception: + error = sys.exc_info()[1] + FreeCAD.Console.PrintError( + "Unexpected error when creating mesh: {}\n" + .format(error) + ) + + doc.recompute() + return doc diff --git a/src/Mod/Fem/femexamples/manager.py b/src/Mod/Fem/femexamples/manager.py index 3b8f330aa8..a0eda86964 100644 --- a/src/Mod/Fem/femexamples/manager.py +++ b/src/Mod/Fem/femexamples/manager.py @@ -71,6 +71,7 @@ def run_all(): run_example("equation_electrostatics_capacitance_two_balls", run_solver=True) run_example("equation_electrostatics_electricforce_elmer_nongui6", run_solver=True) run_example("equation_flow_elmer_2D", run_solver=True) + run_example("equation_flow_turbulent_elmer_2D", run_solver=True) run_example("equation_flux_elmer", run_solver=True) run_example("equation_magnetodynamics_elmer", run_solver=True) run_example("equation_magnetodynamics_2D_elmer.py", run_solver=True) @@ -108,6 +109,7 @@ def setup_all(): run_example("equation_electrostatics_capacitance_two_balls") run_example("equation_electrostatics_electricforce_elmer_nongui6") run_example("equation_flow_elmer_2D") + run_example("equation_flow_turbulent_elmer_2D") run_example("equation_flux_elmer") run_example("equation_magnetodynamics_elmer") run_example("equation_magnetodynamics_2D_elmer.py")