From c06c5788e24bf4daad80367ab089dd4a8ee58aa6 Mon Sep 17 00:00:00 2001 From: Bernd Hahnebach Date: Mon, 1 Aug 2016 21:57:39 +0100 Subject: [PATCH] FEM: constraint thermomech: add implementation for solver CalculiX and fix FEM unit tests --- src/Mod/Fem/FemInputWriter.py | 9 ++ src/Mod/Fem/FemInputWriterCcx.py | 124 +++++++++++++++++- src/Mod/Fem/FemInputWriterZ88.py | 2 + src/Mod/Fem/FemTools.py | 29 +++- src/Mod/Fem/FemToolsCcx.py | 3 +- src/Mod/Fem/FemToolsZ88.py | 1 + src/Mod/Fem/test_files/ccx/cube_frequency.inp | 9 ++ src/Mod/Fem/test_files/ccx/cube_static.inp | 9 ++ 8 files changed, 181 insertions(+), 5 deletions(-) diff --git a/src/Mod/Fem/FemInputWriter.py b/src/Mod/Fem/FemInputWriter.py index 242843bf1..f67975e8f 100644 --- a/src/Mod/Fem/FemInputWriter.py +++ b/src/Mod/Fem/FemInputWriter.py @@ -47,6 +47,7 @@ class FemInputWriter(): fixed_obj, displacement_obj, contact_obj, planerotation_obj, selfweight_obj, force_obj, pressure_obj, + temperature_obj, heatflux_obj, initialtemperature_obj, beamsection_obj, shellthickness_obj, analysis_type, eigenmode_parameters, dir_name @@ -62,6 +63,9 @@ class FemInputWriter(): self.selfweight_objects = selfweight_obj self.force_objects = force_obj self.pressure_objects = pressure_obj + self.temperature_objects = temperature_obj + self.heatflux_objects = heatflux_obj + self.initialtemperature_objects = initialtemperature_obj self.beamsection_objects = beamsection_obj self.shellthickness_objects = shellthickness_obj self.analysis_type = analysis_type @@ -104,6 +108,11 @@ class FemInputWriter(): for femobj in self.planerotation_objects: # femobj --> dict, FreeCAD document object is femobj['Object'] femobj['Nodes'] = FemMeshTools.get_femnodes_by_references(self.femmesh, femobj['Object'].References) + def get_constraints_temperature_nodes(self): + # get nodes + for femobj in self.temperature_objects: # femobj --> dict, FreeCAD document object is femobj['Object'] + femobj['Nodes'] = FemMeshTools.get_femnodes_by_references(self.femmesh, femobj['Object'].References) + def get_constraints_force_nodeloads(self): # check shape type of reference shape for femobj in self.force_objects: # femobj --> dict, FreeCAD document object is femobj['Object'] diff --git a/src/Mod/Fem/FemInputWriterCcx.py b/src/Mod/Fem/FemInputWriterCcx.py index 198f2d79a..2f6e132f2 100644 --- a/src/Mod/Fem/FemInputWriterCcx.py +++ b/src/Mod/Fem/FemInputWriterCcx.py @@ -42,6 +42,7 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): fixed_obj, displacement_obj, contact_obj, planerotation_obj, selfweight_obj, force_obj, pressure_obj, + temperature_obj, heatflux_obj, initialtemperature_obj, beamsection_obj, shellthickness_obj, analysis_type=None, eigenmode_parameters=None, dir_name=None @@ -52,6 +53,7 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): fixed_obj, displacement_obj, contact_obj, planerotation_obj, selfweight_obj, force_obj, pressure_obj, + temperature_obj, heatflux_obj, initialtemperature_obj, beamsection_obj, shellthickness_obj, analysis_type, eigenmode_parameters, dir_name @@ -75,17 +77,28 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): self.write_node_sets_constraints_planerotation(inpfile) if self.contact_objects: self.write_surfaces_contraints_contact(inpfile) + if self.analysis_type == "thermomech" and self.temperature_objects: + self.write_node_sets_constraints_temperature(inpfile) self.write_materials(inpfile) + if self.analysis_type == "thermomech" and self.initialtemperature_objects: + self.write_constraints_initialtemperature(inpfile) self.write_femelementsets(inpfile) if self.planerotation_objects: self.write_constraints_planerotation(inpfile) if self.contact_objects: self.write_constraints_contact(inpfile) - self.write_step_begin(inpfile) + if self.analysis_type == "thermomech": + self.write_step_begin_thermomech(inpfile) + self.write_analysis_thermomech(inpfile) + else: + self.write_step_begin_static_frequency(inpfile) if self.fixed_objects: self.write_constraints_fixed(inpfile) if self.displacement_objects: self.write_constraints_displacement(inpfile) + if self.analysis_type == "thermomech": + self.write_constraints_temperature(inpfile) + self.write_constraints_heatflux(inpfile) if self.analysis_type is None or self.analysis_type == "static": if self.selfweight_objects: self.write_constraints_selfweight(inpfile) @@ -224,16 +237,33 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): for i in v: f.write("{},S{}\n".format(i[0], i[1])) + def write_node_sets_constraints_temperature(self, f): + # get nodes + self.get_constraints_temperature_nodes() + # write nodes to file + f.write('\n***********************************************************\n') + f.write('** Node sets for temperature constraints\n') + f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name)) + for femobj in self.temperature_objects: # femobj --> dict, FreeCAD document object is femobj['Object'] + f.write('*NSET,NSET=' + femobj['Object'].Name + '\n') + for n in femobj['Nodes']: + f.write(str(n) + ',\n') + def write_materials(self, f): f.write('\n***********************************************************\n') f.write('** Materials\n') f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name)) f.write('** Young\'s modulus unit is MPa = N/mm2\n') f.write('** Density\'s unit is t/mm^3\n') + f.write('** Thermal conductivity unit is kW/mm/K = t*mm/K*s^3\n') + f.write('** Specific Heat unit is kJ/t/K = mm^2/s^2/K\n') for femobj in self.material_objects: # femobj --> dict, FreeCAD document object is femobj['Object'] mat_obj = femobj['Object'] # get material properties - Currently in SI units: M/kg/s/Kelvin YM_in_MPa = 1 + TC_in_WmK = 1 + TEC_in_mmK = 1 + SH_in_JkgK = 1 PR = 1 density_in_tonne_per_mm3 = 1 try: @@ -245,6 +275,21 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): PR = float(mat_obj.Material['PoissonRatio']) except: FreeCAD.Console.PrintError("No PoissonRatio defined for material: default used\n") + try: + TC = FreeCAD.Units.Quantity(mat_obj.Material['ThermalConductivity']) + TC_in_WmK = TC.getValueAs('W/m/K') # SvdW: Add factor to force units to results' base units of t/mm/s/K - W/m/K results in no factor needed + except: + FreeCAD.Console.PrintError("No ThermalConductivity defined for material: default used\n") + try: + TEC = FreeCAD.Units.Quantity(mat_obj.Material['ThermalExpansionCoefficient']) + TEC_in_mmK = TEC.getValueAs('mm/mm/K') + except: + FreeCAD.Console.PrintError("No ThermalExpansionCoefficient defined for material: default used\n") + try: + SH = FreeCAD.Units.Quantity(mat_obj.Material['SpecificHeat']) + SH_in_JkgK = SH.getValueAs('J/kg/K') * 1e+06 # SvdW: Add factor to force units to results' base units of t/mm/s/K + except: + FreeCAD.Console.PrintError("No SpecificHeat defined for material: default used\n") mat_info_name = mat_obj.Material['Name'] mat_name = mat_obj.Name # write material properties @@ -259,6 +304,12 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): FreeCAD.Console.PrintError("No Density defined for material: default used\n") f.write('*DENSITY \n') f.write('{0:.3e}, \n'.format(density_in_tonne_per_mm3)) + f.write('*CONDUCTIVITY \n') + f.write('{}, \n'.format(TC_in_WmK)) + f.write('*EXPANSION \n') + f.write('{}, \n'.format(TEC_in_mmK)) + f.write('*SPECIFIC HEAT \n') + f.write('{}, \n'.format(SH_in_JkgK)) def write_femelementsets(self, f): f.write('\n***********************************************************\n') @@ -295,7 +346,7 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): setion_def = '*SOLID SECTION, ' + elsetdef + material + '\n' f.write(setion_def) - def write_step_begin(self, f): + def write_step_begin_static_frequency(self, f): f.write('\n***********************************************************\n') f.write('** One step is needed to run the mechanical analysis of FreeCAD\n') f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name)) @@ -316,6 +367,15 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): analysis_static += ', SOLVER=ITERATIVE CHOLESKY' f.write(analysis_static + '\n') + def write_step_begin_thermomech(self, f): + f.write('\n***********************************************************\n') + f.write('** One step is needed to run the thermomechanical analysis of FreeCAD\n') + f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name)) + thermomech_step = '*STEP' + if self.solver_obj.GeometricalNonlinearity == "nonlinear": + thermomech_step += ', NLGEOM' + thermomech_step += ', INC=' + str(self.solver_obj.IterationsMaximum) + f.write(thermomech_step + '\n') def write_constraints_fixed(self, f): f.write('\n***********************************************************\n') @@ -399,6 +459,15 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): fric_obj_name = femobj['Object'].Name f.write('*MPC\n') f.write('PLANE,' + fric_obj_name + '\n') + + def write_constraints_temperature(self, f): + f.write('\n***********************************************************\n') + f.write('** Fixed temperature constraint applied\n') + f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name)) + for ftobj in self.temperature_objects: + fixedtemp_obj = ftobj['Object'] + f.write('*BOUNDARY\n') + f.write('{},11,11,{}\n'.format(fixedtemp_obj.Name, fixedtemp_obj.Temperature)) f.write('\n') def write_constraints_selfweight(self, f): @@ -459,6 +528,22 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): for i in v: f.write("{},P{},{}\n".format(i[0], i[1], rev * prs_obj.Pressure)) + def write_constraints_heatflux(self, f): + f.write('\n***********************************************************\n') + f.write('** Heatflux constraints\n') + f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name)) + for hfobj in self.heatflux_objects: + heatflux_obj = hfobj['Object'] + f.write('*FILM\n') + for o, elem_tup in heatflux_obj.References: + for elem in elem_tup: + ho = o.Shape.getElement(elem) + if ho.ShapeType == 'Face': + v = self.mesh_object.FemMesh.getccxVolumesByFace(ho) + f.write("** Heat flux on face {}\n".format(elem)) + for i in v: + f.write("{},F{},{},{}\n".format(i[0], i[1], heatflux_obj.AmbientTemp, heatflux_obj.FilmCoef * 0.001)) # SvdW add factor to force heatflux to units system of t/mm/s/K # OvG: Only write out the VolumeIDs linked to a particular face + def write_analysis_frequency(self, f): f.write('\n***********************************************************\n') f.write('** Frequency analysis\n') @@ -466,6 +551,36 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): f.write('*FREQUENCY\n') f.write('{},{},{}\n'.format(self.no_of_eigenfrequencies, self.eigenfrequeny_range_low, self.eigenfrequeny_range_high)) + def write_analysis_thermomech(self, f): + f.write('\n***********************************************************\n') + f.write('** Un-Coupled temperature displacement analysis\n') + f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name)) + thermomech_analysis = '*COUPLED TEMPERATURE-DISPLACEMENT' + if self.solver_obj.MatrixSolverType == "default": + pass + elif self.solver_obj.MatrixSolverType == "spooles": + thermomech_analysis += ', SOLVER=SPOOLES' + elif self.solver_obj.MatrixSolverType == "iterativescaling": + thermomech_analysis += ', SOLVER=ITERATIVE SCALING' + elif self.solver_obj.MatrixSolverType == "iterativecholesky": + thermomech_analysis += ', SOLVER=ITERATIVE CHOLESKY' + if self.solver_obj.SteadyState: + thermomech_analysis += ', STEADY STATE\n' + self.solver_obj.TimeInitialStep = 1.0 # Set time to 1 and ignore user inputs for steady state + self.solver_obj.TimeEnd = 1.0 + thermomech_time = '{},{}'.format(self.solver_obj.TimeInitialStep, self.solver_obj.TimeEnd) # OvG: 1.0 increment, total time 1 for steady state will cut back automatically + f.write(thermomech_analysis + '\n') + f.write(thermomech_time + '\n') + + def write_constraints_initialtemperature(self, f): + f.write('\n***********************************************************\n') + f.write('** Initial temperature constraint\n') + f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name)) + f.write('*INITIAL CONDITIONS,TYPE=TEMPERATURE\n') + for itobj in self.initialtemperature_objects: # Should only be one + inittemp_obj = itobj['Object'] + f.write('Nall,{}\n'.format(inittemp_obj.initialTemperature)) # OvG: Initial temperature + def write_outputs_types(self, f): f.write('\n***********************************************************\n') f.write('** Outputs --> frd file\n') @@ -474,7 +589,10 @@ class FemInputWriterCcx(FemInputWriter.FemInputWriter): f.write('*NODE FILE, OUTPUT=2d\n') else: f.write('*NODE FILE\n') - f.write('U\n') + if self.analysis_type == "thermomech": # MPH write out nodal temperatures if thermomechanical + f.write('U, NT\n') + else: + f.write('U\n') f.write('*EL FILE\n') f.write('S, E\n') f.write('** outputs --> dat file\n') diff --git a/src/Mod/Fem/FemInputWriterZ88.py b/src/Mod/Fem/FemInputWriterZ88.py index 60a7abfaf..ccf4af1df 100644 --- a/src/Mod/Fem/FemInputWriterZ88.py +++ b/src/Mod/Fem/FemInputWriterZ88.py @@ -39,6 +39,7 @@ class FemInputWriterZ88(FemInputWriter.FemInputWriter): fixed_obj, displacement_obj, contact_obj, planerotation_obj, selfweight_obj, force_obj, pressure_obj, + temperature_obj, heatflux_obj, initialtemperature_obj, beamsection_obj, shellthickness_obj, analysis_type=None, eigenmode_parameters=None, dir_name=None @@ -49,6 +50,7 @@ class FemInputWriterZ88(FemInputWriter.FemInputWriter): fixed_obj, displacement_obj, contact_obj, planerotation_obj, selfweight_obj, force_obj, pressure_obj, + temperature_obj, heatflux_obj, initialtemperature_obj, beamsection_obj, shellthickness_obj, analysis_type, eigenmode_parameters, dir_name diff --git a/src/Mod/Fem/FemTools.py b/src/Mod/Fem/FemTools.py index e6c90967b..38748999d 100644 --- a/src/Mod/Fem/FemTools.py +++ b/src/Mod/Fem/FemTools.py @@ -151,6 +151,9 @@ class FemTools(QtCore.QRunnable, QtCore.QObject): # [{'Object':fixed_constraints, 'NodeSupports':bool}, {}, ...] # [{'Object':force_constraints, 'NodeLoad':value}, {}, ... # [{'Object':pressure_constraints, 'xxxxxxxx':value}, {}, ...] + # [{'Object':temerature_constraints, 'xxxxxxxx':value}, {}, ...] + # [{'Object':heatflux_constraints, 'xxxxxxxx':value}, {}, ...] + # [{'Object':initialtemperature_constraints, 'xxxxxxxx':value}, {}, ...] # [{'Object':beam_sections, 'xxxxxxxx':value}, {}, ...] # [{'Object':shell_thicknesses, 'xxxxxxxx':value}, {}, ...] # [{'Object':contact_constraints, 'xxxxxxxx':value}, {}, ...] @@ -190,6 +193,18 @@ class FemTools(QtCore.QRunnable, QtCore.QObject): # set of displacements for the analysis. Updated with update_objects # Individual displacement_constraints are Proxy.Type "FemConstraintDisplacement" self.displacement_constraints = [] + ## @var temperature_constraints + # set of temperatures for the analysis. Updated with update_objects + # Individual temperature_constraints are Proxy.Type "FemConstraintTemperature" + self.temperature_constraints = [] + ## @var heatflux_constraints + # set of heatflux constraints for the analysis. Updated with update_objects + # Individual heatflux_constraints are Proxy.Type "FemConstraintHeatflux" + self.heatflux_constraints = [] + ## @var initialtemperature_constraints + # set of initial temperatures for the analysis. Updated with update_objects + # Individual initialTemperature_constraints are Proxy.Type "FemConstraintInitialTemperature" + self.initialtemperature_constraints = [] ## @var planerotation_constraints # set of plane rotation constraints from the analysis. Updated with update_objects # Individual constraints are "Fem::ConstraintPlaneRotation" type @@ -240,10 +255,22 @@ class FemTools(QtCore.QRunnable, QtCore.QObject): PressureObjectDict = {} PressureObjectDict['Object'] = m self.pressure_constraints.append(PressureObjectDict) - elif m.isDerivedFrom("Fem::ConstraintDisplacement"): # OvG: Replacement reference to C++ implementation of Displacement Constraint + elif m.isDerivedFrom("Fem::ConstraintDisplacement"): displacement_constraint_dict = {} displacement_constraint_dict['Object'] = m self.displacement_constraints.append(displacement_constraint_dict) + elif m.isDerivedFrom("Fem::ConstraintTemperature"): + temperature_constraint_dict = {} + temperature_constraint_dict['Object'] = m + self.temperature_constraints.append(temperature_constraint_dict) + elif m.isDerivedFrom("Fem::ConstraintHeatflux"): + heatflux_constraint_dict = {} + heatflux_constraint_dict['Object'] = m + self.heatflux_constraints.append(heatflux_constraint_dict) + elif m.isDerivedFrom("Fem::ConstraintInitialTemperature"): + initialtemperature_constraint_dict = {} + initialtemperature_constraint_dict['Object'] = m + self.initialtemperature_constraints.append(initialtemperature_constraint_dict) elif m.isDerivedFrom("Fem::ConstraintPlaneRotation"): planerotation_constraint_dict = {} planerotation_constraint_dict['Object'] = m diff --git a/src/Mod/Fem/FemToolsCcx.py b/src/Mod/Fem/FemToolsCcx.py index 56f374dd1..517865aaa 100644 --- a/src/Mod/Fem/FemToolsCcx.py +++ b/src/Mod/Fem/FemToolsCcx.py @@ -34,7 +34,7 @@ from PySide import QtCore class FemToolsCcx(FemTools.FemTools): - known_analysis_types = ["static", "frequency"] + known_analysis_types = ["static", "frequency", "thermomech"] finished = QtCore.Signal(int) ## The constructor @@ -93,6 +93,7 @@ class FemToolsCcx(FemTools.FemTools): self.fixed_constraints, self.displacement_constraints, self.contact_constraints, self.planerotation_constraints, self.selfweight_constraints, self.force_constraints, self.pressure_constraints, + self.temperature_constraints, self.heatflux_constraints, self.initialtemperature_constraints, self.beam_sections, self.shell_thicknesses, self.analysis_type, self.eigenmode_parameters, self.working_dir diff --git a/src/Mod/Fem/FemToolsZ88.py b/src/Mod/Fem/FemToolsZ88.py index 083a28bb6..bc965bbac 100644 --- a/src/Mod/Fem/FemToolsZ88.py +++ b/src/Mod/Fem/FemToolsZ88.py @@ -86,6 +86,7 @@ class FemToolsZ88(FemTools.FemTools): self.fixed_constraints, self.displacement_constraints, self.contact_constraints, self.planerotation_constraints, self.selfweight_constraints, self.force_constraints, self.pressure_constraints, + self.temperature_constraints, self.heatflux_constraints, self.initialtemperature_constraints, self.beam_sections, self.shell_thicknesses, self.analysis_type, None, self.working_dir diff --git a/src/Mod/Fem/test_files/ccx/cube_frequency.inp b/src/Mod/Fem/test_files/ccx/cube_frequency.inp index e3ee6565f..c2dbe8d41 100644 --- a/src/Mod/Fem/test_files/ccx/cube_frequency.inp +++ b/src/Mod/Fem/test_files/ccx/cube_frequency.inp @@ -468,12 +468,21 @@ Eall ** Materials ** written by write_materials function ** Young's modulus unit is MPa = N/mm2 +** Density's unit is t/mm^3 +** Thermal conductivity unit is kW/mm/K = t*mm/K*s^3 +** Specific Heat unit is kJ/t/K = mm^2/s^2/K **FreeCAD material name: Steel-Generic *MATERIAL, NAME=MechanicalMaterial *ELASTIC 200000 , 0.300 *DENSITY 7.900e-09, +*CONDUCTIVITY +1, +*EXPANSION +1, +*SPECIFIC HEAT +1, *********************************************************** ** Sections diff --git a/src/Mod/Fem/test_files/ccx/cube_static.inp b/src/Mod/Fem/test_files/ccx/cube_static.inp index e8edf9e0d..14765d2e3 100644 --- a/src/Mod/Fem/test_files/ccx/cube_static.inp +++ b/src/Mod/Fem/test_files/ccx/cube_static.inp @@ -468,12 +468,21 @@ Eall ** Materials ** written by write_materials function ** Young's modulus unit is MPa = N/mm2 +** Density's unit is t/mm^3 +** Thermal conductivity unit is kW/mm/K = t*mm/K*s^3 +** Specific Heat unit is kJ/t/K = mm^2/s^2/K **FreeCAD material name: Steel-Generic *MATERIAL, NAME=MechanicalMaterial *ELASTIC 200000 , 0.300 *DENSITY 7.900e-09, +*CONDUCTIVITY +1, +*EXPANSION +1, +*SPECIFIC HEAT +1, *********************************************************** ** Sections