Scripting setup of boundary conditions

[markrmau]

Firstly thanks for donig such a great job with this tool.

I used scripting in freecad to setup a relatively complex geometry with many faces. Ideally I'd like to be able to name the faces during the model generation (for example "walls") and then set these as noslip BCs for the OF simulation. However, I see this is not possible yet (https://www.freecadweb.org/wiki/Naming_project)

Another alternative would be to run a script that reads the mesh files, determines which are the faces I require for walls and inlet / outlet and writes it into the freecad file CfdAnalysis object - but this scripting doesn't seem to be possible either?

A 3rd alternative is to write a script which re-writes createPatchDict before runing the simulation. Is this the way it needs to be done?

Thanks for any comments.

[chrisb]

Hi and welcome to the forum!

You may have a look at realthunders branch. He has developed techniques to master (at least parts of) the topological naming issues. That may be helpful.

[herbk]

Hi Mark,

If you have many faces, the most will be "no slip" i think, use the option "select from list" of the UI. If using it, use the "select all" option and deselect inlet, outlet and walls with slip (if you have), which are usualy a few.

My order to select: inlet, outlet and walls with slip, so i have the faces nubers for it.

[bernd]

Another alternative would be to run a script that reads the mesh files, determines which are the faces I require for walls and inlet / outlet and writes it into the freecad file CfdAnalysis object - but this scripting doesn't seem to be possible either?

this is possible for FEM. Since Cfd is based on former Fem code it should be possible for Cfd either. If you would post a specific (but simple) example we may could set this up together.

[markrmau]

Thanks for the many prompt replies!

I will try to compile realthinders branch at some point later as I'm sure that's a better way to go.

My problem is I have approx 60 faces to be added to Wall BC and about 60 to be added to the Outlet BC

I think I see how to script the addition now - thanks bernd. I didn't realise python console would show me the script as executed.

The script below will setup the part, add the cfd analysis and add a BC of Face138 to wall. So that should get me going. Note I'm sure more clever coding could be done to create the object, but this works.

I now need to work out a way to find the face names of the outer circles on the dimples. One way would be to write a mesh geometry, and reverse engineer the face locations but I'll try for something better!

Thanks for the help.

Code: Select all

import Part, math
from FreeCAD import Base

# t-piece radius and height
r=75.0
H=1000.0

# dimple separation, height and radius
h_sep = 100
d_h=6.0
d_r=2.5

mydoc = FreeCAD.newDocument("t-piece")
cyl_in_H = Part.makeCylinder(r,H/2,Base.Vector(-H/2,0,0),Base.Vector(1,0,0),360)
cyl_in_V = Part.makeCylinder(r,H,Base.Vector(0,0,-H/2),Base.Vector(0,0,1),360)
inlet = cyl_in_H.fuse(cyl_in_V)

for th in range(0,360,30):
	x=(r-d_h/2)*math.cos(th*math.pi/180)
	y=(r-d_h/2)*math.sin(th*math.pi/180)

	h=r+h_sep
	while h <= H/2-h_sep:	
		dimple = Part.makeCylinder(d_r,d_h,Base.Vector(x,y,h),Base.Vector(x,y,0),360)
		inlet=inlet.fuse(dimple)
		dimple = Part.makeCylinder(d_r,d_h,Base.Vector(x,y,-h),Base.Vector(x,y,0),360)
		inlet=inlet.fuse(dimple)
		h=h+h_sep

Part.show(inlet)
mydoc.recompute()
Gui.SendMsgToActiveView("ViewFit")

### Begin command Cfd_Analysis
analysis = CfdAnalysis.makeCfdAnalysis('CfdAnalysis')
CfdTools.setActiveAnalysis(analysis)
analysis.addObject(CfdPhysicsSelection.makeCfdPhysicsSelection())
analysis.addObject(CfdFluidMaterial.makeCfdFluidMaterial('FluidProperties'))
analysis.addObject(CfdInitialiseFlowField.makeCfdInitialFlowField())
analysis.addObject(CfdSolverFoam.makeCfdSolverFoam())
### End command Cfd_Analysis
# Gui.Selection.addSelection('t_piece','CfdAnalysis')
App.activeDocument().recompute(None,True,True)
### Begin command Cfd_FluidBoundary

CfdTools.getActiveAnalysis().addObject(CfdFluidBoundary.makeCfdFluidBoundary())
### End command Cfd_FluidBoundary
# Gui.Selection.clearSelection()

bc = FreeCAD.ActiveDocument.CfdFluidBoundary
bc.BoundaryType = 'wall'
bc.BoundarySubType = 'fixedWall'
bc.ThermalBoundaryType = 'fixedValue'
bc.VelocityIsCartesian = True
bc.Ux = '0 mm/s'
bc.Uy = '0 mm/s'
bc.Uz = '0 mm/s'
bc.VelocityMag = '0 mm/s'
bc.DirectionFace = ''
bc.ReverseNormal = False
bc.MassFlowRate = '0 kg/s'
bc.VolFlowRate = '0 mm^3/s'
bc.Pressure = '0 kg/(mm*s^2)'
bc.SlipRatio = '0 '
bc.Temperature = '290 K'
bc.HeatFlux = '0 kg/s^3'
bc.HeatTransferCoeff = '0 kg/(s^3*K)'
bc.TurbulenceInletSpecification = 'intensityAndLengthScale'
bc.TurbulentKineticEnergy = '10000 mm^2/s^2'
bc.SpecificDissipationRate = '57 deg/s'
bc.TurbulenceIntensity = '0.1 '
bc.TurbulenceLengthScale = '100 mm'
bc.VolumeFractions = {}
bc.PorousBaffleMethod = 'porousCoeff'
bc.PressureDropCoeff = '0 '
bc.ScreenWireDiameter = '0.2 mm'
bc.ScreenSpacing = '2 mm'
FreeCAD.ActiveDocument.CfdFluidBoundary.Label = 'wall'
FreeCAD.ActiveDocument.CfdFluidBoundary.References = []
FreeCAD.ActiveDocument.CfdFluidBoundary.References.append(('Shape', 'Face138'))
FreeCAD.ActiveDocument.recompute()

[markrmau]

Oh, I understand now - silly me!

'Face138' (my inlet) is the entry number 138 (if counted from 1) in the Faces list.

So

Code: Select all

f=App.ActiveDocument.Shape.Shape.Faces[137]
f.CenterOfMass

will give me the center of the inlet face (within a tolerance).

So after the shape is recomputed the BC's can be assigned by iterating the Faces list.

[bernd]

Yep you got it ...

this may help with improving your script FEM_Tutorial_Python But keep in mind cfd is based on former Fem code, means something may work slightly different.

[markrmau]

Thanks for the tutorial !

I have got it to the point where I want it - so that everything is setup ready to write and check the mesh and solver case files - a point where review is probably needed before charging ahead.

There may be errors in selection of faces still - I haven't fully checked. Perhaps I should tidy this example up (break into functions) and make as a scripting tutorial for cfdof? (comments welcome).

I haven't got a good working mesh and solution yet - will work on it later.

My intention is to use a gmesh. I suspect the inlet BC I have put in the macro (volumetricFlowRateInlet) has an error (in the openfoam code). I'll check this though - it might be my error.

Thanks again, Mark.

Code: Select all

import Part, math
import CfdAnalysis, CfdTools, CfdPhysicsSelection,CfdFluidMaterial, CfdInitialiseFlowField, CfdSolverFoam, CfdFluidBoundary, CfdMesh
from FreeCAD import Base

# t-piece radius and height
r=75.0
H=1000.0

# dimple separation, height and radius
h_sep = 50
d_h=8.0
d_r=5

mydoc = FreeCAD.newDocument('t-piece')
tp = mydoc.addObject('Part::Feature','tp')
cyl_in_H = Part.makeCylinder(r,H/2,Base.Vector(-H/2,0,0),Base.Vector(1,0,0),360)
cyl_in_V = Part.makeCylinder(r,H,Base.Vector(0,0,-H/2),Base.Vector(0,0,1),360)
tp.Shape = cyl_in_H.fuse(cyl_in_V)

for th in range(0,360,15):
	x=(r-d_h/2)*math.cos(th*math.pi/180)
	y=(r-d_h/2)*math.sin(th*math.pi/180)

	h=r+h_sep
	while h <= H/2-h_sep:	
		dimple = Part.makeCylinder(d_r,d_h,Base.Vector(x,y,h),Base.Vector(x,y,0),360)
		tp.Shape=tp.Shape.fuse(dimple)
		dimple = Part.makeCylinder(d_r,d_h,Base.Vector(x,y,-h),Base.Vector(x,y,0),360)
		tp.Shape=tp.Shape.fuse(dimple)
		h=h+h_sep

mydoc.recompute()

# View part
import FreeCADGui
FreeCADGui.ActiveDocument.activeView().viewAxonometric()
FreeCADGui.SendMsgToActiveView('ViewFit')

# Create analysis container
analysis_object = CfdAnalysis.makeCfdAnalysis('CfdAnalysis')

# Setup cfd physics model
cfd_phys_sel = CfdPhysicsSelection.makeCfdPhysicsSelection()
cfd_phys_sel.Time = 'Steady'
cfd_phys_sel.Phase = 'Single'
cfd_phys_sel.Flow = 'Incompressible'
cfd_phys_sel.Thermal = 'None'
cfd_phys_sel.Turbulence = 'RANS'
cfd_phys_sel.TurbulenceModel = 'kOmegaSST'
cfd_phys_sel.gx = '0 mm/s^2'
cfd_phys_sel.gy = '-9.8e+03 mm/s^2'
cfd_phys_sel.gz = '0 mm/s^2'

# Setup fluid (water)
cfd_fluid = CfdFluidMaterial.makeCfdFluidMaterial('FluidProperties')
cfd_fluid.Material = {'CardName': 'Water', 'AuthorAndLicense': 'Water', 'Name': 'Water', \
	'Description': 'Standard distilled water properties at 20 Degrees Celsius and 1 atm', \
	'Density': '998 kg/m^3', 'DynamicViscosity': '1.003e-3 kg/m/s'}

# Setup initialisation (leave defaults)
cfd_init = CfdInitialiseFlowField.makeCfdInitialFlowField()

# Create solver object
cfd_solver = CfdSolverFoam.makeCfdSolverFoam()

# Setup boundary conditions.  Note all BC options are included but only those relevant to the BC subtype are used
wall_BC = CfdFluidBoundary.makeCfdFluidBoundary()
wall_BC.BoundaryType = 'wall'
wall_BC.BoundarySubType = 'fixedWall'
wall_BC.ThermalBoundaryType = 'fixedValue'
wall_BC.VelocityIsCartesian = True
wall_BC.Ux = '0 mm/s'
wall_BC.Uy = '0 mm/s'
wall_BC.Uz = '0 mm/s'
wall_BC.VelocityMag = '0 mm/s'
wall_BC.DirectionFace = ''
wall_BC.ReverseNormal = False
wall_BC.MassFlowRate = '0 kg/s'
wall_BC.VolFlowRate = '0 mm^3/s'
wall_BC.Pressure = '0 kg/(mm*s^2)'
wall_BC.SlipRatio = '0 '
wall_BC.Temperature = '290 K'
wall_BC.HeatFlux = '0 kg/s^3'
wall_BC.HeatTransferCoeff = '0 kg/(s^3*K)'
wall_BC.TurbulenceInletSpecification = 'intensityAndLengthScale'
wall_BC.TurbulentKineticEnergy = '10000 mm^2/s^2'
wall_BC.SpecificDissipationRate = '57 deg/s'
wall_BC.TurbulenceIntensity = '0.1 '
wall_BC.TurbulenceLengthScale = '100 mm'
wall_BC.VolumeFractions = {}
wall_BC.PorousBaffleMethod = 'porousCoeff'
wall_BC.PressureDropCoeff = '0 '
wall_BC.ScreenWireDiameter = '0.2 mm'
wall_BC.ScreenSpacing = '2 mm'
wall_BC.References = []
wall_BC.Label = 'wall'

inlet_BC = CfdFluidBoundary.makeCfdFluidBoundary()
inlet_BC.BoundaryType = 'inlet'
inlet_BC.BoundarySubType = 'volumetricFlowRateInlet'
inlet_BC.ThermalBoundaryType = 'fixedValue'
inlet_BC.VelocityIsCartesian = True
inlet_BC.Ux = '0 mm/s'
inlet_BC.Uy = '0 mm/s'
inlet_BC.Uz = '0 mm/s'
inlet_BC.VelocityMag = '0 mm/s'
inlet_BC.DirectionFace = ''
inlet_BC.ReverseNormal = True
inlet_BC.MassFlowRate = '0 kg/s'
inlet_BC.VolFlowRate = '1000 mm^3/s'
inlet_BC.Pressure = '0 kg/(mm*s^2)'
inlet_BC.SlipRatio = '0 '
inlet_BC.Temperature = '290 K'
inlet_BC.HeatFlux = '0 kg/s^3'
inlet_BC.HeatTransferCoeff = '0 kg/(s^3*K)'
inlet_BC.TurbulenceInletSpecification = 'intensityAndLengthScale'
inlet_BC.TurbulentKineticEnergy = '10000 mm^2/s^2'
inlet_BC.SpecificDissipationRate = '57 deg/s'
inlet_BC.TurbulenceIntensity = '0.1 '
inlet_BC.TurbulenceLengthScale = '100 mm'
inlet_BC.VolumeFractions = {}
inlet_BC.PorousBaffleMethod = 'porousCoeff'
inlet_BC.PressureDropCoeff = '0 '
inlet_BC.ScreenWireDiameter = '0.2 mm'
inlet_BC.ScreenSpacing = '2 mm'
inlet_BC.References = []
inlet_BC.Label = 'inlet'

outlet_BC = CfdFluidBoundary.makeCfdFluidBoundary()
outlet_BC.BoundaryType = 'outlet'
outlet_BC.BoundarySubType = 'staticPressureOutlet'
outlet_BC.ThermalBoundaryType = 'fixedValue'
outlet_BC.VelocityIsCartesian = True
outlet_BC.Ux = '0 mm/s'
outlet_BC.Uy = '0 mm/s'
outlet_BC.Uz = '0 mm/s'
outlet_BC.VelocityMag = '0 mm/s'
outlet_BC.DirectionFace = ''
outlet_BC.ReverseNormal = False
outlet_BC.MassFlowRate = '0 kg/s'
outlet_BC.VolFlowRate = '0 mm^3/s'
outlet_BC.Pressure = '0 kg/(mm*s^2)'
outlet_BC.SlipRatio = '0 '
outlet_BC.Temperature = '290 K'
outlet_BC.HeatFlux = '0 kg/s^3'
outlet_BC.HeatTransferCoeff = '0 kg/(s^3*K)'
outlet_BC.TurbulenceInletSpecification = 'intensityAndLengthScale'
outlet_BC.TurbulentKineticEnergy = '10000 mm^2/s^2'
outlet_BC.SpecificDissipationRate = '57 deg/s'
outlet_BC.TurbulenceIntensity = '0.1 '
outlet_BC.TurbulenceLengthScale = '100 mm'
outlet_BC.VolumeFractions = {}
outlet_BC.PorousBaffleMethod = 'porousCoeff'
outlet_BC.PressureDropCoeff = '0 '
outlet_BC.ScreenWireDiameter = '0.2 mm'
outlet_BC.ScreenSpacing = '2 mm'
outlet_BC.References = []
outlet_BC.Label = 'outlet'

# Traverse faces list and assign faces to appropriate BC
tol = 1e-6
for i,tp_face in enumerate(tp.Shape.Faces):
	com = tp_face.CenterOfMass
	if (com-Base.Vector(-H/2,0,0)).Length<tol:
		# The inlet is at x=-H/2
		inlet_BC.References.append(('tp','Face'+str(i+1)))
		print('Inlet: '+ 'Face'+str(i+1))
	elif len(tp_face.Edges)<2 and abs(com[2])<H/2-tol and (math.sqrt(com[0]**2+com[1]**2)>r+d_h/2-tol):
		# The pipe from inlet to t has 3 face, and 
		# the outlets are not the ends of the t, 
		# but are at a radius of r+d_h/2
		outlet_BC.References.append(('tp','Face'+str(i+1)))
		print('Outlet: '+ 'Face'+str(i+1))
		
	else:
		# Otherwise it is a wall
		wall_BC.References.append(('tp','Face'+str(i+1)))
		print('Wall: '+ 'Face'+str(i+1))

# Add objects to analysis container
analysis_object.addObject(cfd_phys_sel)
analysis_object.addObject(cfd_fluid)
analysis_object.addObject(cfd_init)
analysis_object.addObject(cfd_solver)
analysis_object.addObject(wall_BC)
analysis_object.addObject(inlet_BC)
analysis_object.addObject(outlet_BC)

CfdTools.setActiveAnalysis(analysis_object)

tp_mesh = CfdMesh.makeCfdMesh('tp_Mesh')
tp_mesh.Part=tp
analysis_object.addObject(tp_mesh)

App.activeDocument().recompute(None,True,True)

[bernd]

Perhaps I should tidy this example up (break into functions) and make as a scripting tutorial for cfdof? (comments welcome).

a scripting tutorial for cfd would be cool. IMO functions are not needed. This makes it more difficault to understan for beginners. If there are no functions the code needs just to be copied in to Python console.

[oliveroxtoby]

There may be errors in selection of faces still - I haven't fully checked. Perhaps I should tidy this example up (break into functions) and make as a scripting tutorial for cfdof? (comments welcome).

Any tutorial/example/training material would be a great contribution.