Help Wanted - Arbitrary (general) beam section in FEM workbench

About the development of the FEM module/workbench.

Moderator: bernd

Post Reply
User avatar
NewJoker
Veteran
Posts: 3016
Joined: Sun Oct 11, 2020 7:49 pm

Help Wanted - Arbitrary (general) beam section in FEM workbench

Post by NewJoker »

Hi,

here we recently had a short discussion about the general beam section (CalculiX features) that allows users to define beam elements with arbitrary cross-sections: https://forum.freecadweb.org/viewtopic. ... 4&start=20

In my opinion, it could be very useful so I'd like to continue this topic and move on to the potential implementation.

Let's start from what is required by CalculiX (apart from standard definitions) to use this feature (source: http://www.dhondt.de/ccx_2.17.pdf):

Code: Select all

*ELEMENT, TYPE=U1, ...
...
*USER ELEMENT, TYPE=U1, NODES=2, INTEGRATION POINTS=2, MAXDOF=6
*BEAM SECTION,SECTION=GENERAL, ...
A, I_11, I_12, I_22, k
uv_x, uv_y, uv_z
where: A - cross-section area, I_11, I_12 and I_22 - moments of inertia (I_12 must be zero), k - Timoshenko shear coefficient, uv_x, uv_y and uv_z - coordinates of a unit vector in 1 direction (they are used for other section types as well but for this one they are required).

Now let's get to my implementation attempt. Please keep in mind that it's my first time trying to add code to FreeCAD so I may make some mistakes.
From what I've found, at least 3 files have to be edited to add new section type. Below are the links to those files in github repository together with my suggestion of what to add to them:

write_femelement_geometry.py
https://github.com/FreeCAD/FreeCAD/blob ... eometry.py

Code: Select all

if beamsec_obj.SectionType == "General":
                    area = beamsec_obj.GenArea.getValueAs("mm^2").Value
                    in11 = beamsec_obj.GenIn11.getValueAs("mm^4").Value
                    in12 = 0
                    in22 = beamsec_obj.GenIn22.getValueAs("mm^4").Value
                    k = beamsec_obj.GenK.getValueAs("mm^4").Value
                    section_type = ", SECTION=GEN"
                    section_geo = "{:.13G},{:.13G},{:.13G},{:.13G},{:.13G}\n".format(area, in11, 0 , in22, k)
                    section_def = "*BEAM SECTION, {}{}{}\n".format(
                        elsetdef,
                        material,
                        section_type
                    )

element_geometry1D.py
https://github.com/FreeCAD/FreeCAD/blob ... metry1D.py

Code: Select all

 Type = "Fem::ElementGeometry1D"
    known_beam_types = ["Rectangular", "Circular", "Pipe", "General"]

    def __init__(self, obj):
        super(ElementGeometry1D, self).__init__(obj)

        obj.addProperty(
            "App::PropertyLength",
            "GenArea",
            "GenBeamSection",
            "set area of the general beam elements"
        )

        obj.addProperty(
            "App::PropertyLength",
            "GenIn11",
            "GenBeamSection",
            "set I_11 of the general beam elements"
        )
        
          obj.addProperty(
            "App::PropertyLength",
            "GenIn22",
            "GenBeamSection",
            "set I_22 of the general beam elements"
        )
        
          obj.addProperty(
            "App::PropertyLength",
            "GenK",
            "GenBeamSection",
            "set Timoshenko shear coefficient of the general beam elements"
        )
task_element_geometry1D.py
https://github.com/FreeCAD/FreeCAD/blob ... metry1D.py

Code: Select all

QtCore.QObject.connect(
            self.parameterWidget.if_gen_area,
            QtCore.SIGNAL("valueChanged(Base::Quantity)"),
            self.gen_area_changed
        )
            QtCore.QObject.connect(
            self.parameterWidget.if_gen_in11,
            QtCore.SIGNAL("valueChanged(Base::Quantity)"),
            self.gen_in11_changed
        )
        QtCore.QObject.connect(
            self.parameterWidget.if_gen_in22,
            QtCore.SIGNAL("valueChanged(Base::Quantity)"),
            self.gen_in22_changed
        )
        QtCore.QObject.connect(
            self.parameterWidget.if_gen_k,
            QtCore.SIGNAL("valueChanged(Base::Quantity)"),
            self.gen_k_changed
        )
...
...
...        
   def get_beamsection_props(self):
        self.GenArea = self.obj.GenArea
        self.GenIn11 = self.obj.GenIn11
        self.GenIn22 = self.obj.GenIn22
        self.GenK = self.obj.GenK
   def set_beamsection_props(self):
        self.GenArea = self.obj.GenArea
        self.GenIn11 = self.obj.GenIn11
        self.GenIn22 = self.obj.GenIn22
        self.GenK = self.obj.GenK
   def updateParameterWidget(self):
        "fills the widgets"
        index_crosssectiontype = self.parameterWidget.cb_crosssectiontype.findText(
            self.SectionType
        )
        self.parameterWidget.cb_crosssectiontype.setCurrentIndex(index_crosssectiontype)
        self.parameterWidget.if_gen_area.setText(self.GenArea.UserString)
        self.parameterWidget.if_gen_in11.setText(self.GenIn11.UserString)
        self.parameterWidget.if_gen_in22.setText(self.GenIn22.UserString)
        self.parameterWidget.if_gen_k.setText(self.GenK.UserString)
...
...
...           
   def gen_area_changed(self, base_quantity_value):
        self.GenArea = base_quantity_value
       
   def gen_in11_changed(self, base_quantity_value):
        self.GenIn11 = base_quantity_value
        
   def gen_in22_changed(self, base_quantity_value):
        self.GenIn22 = base_quantity_value
        
   def gen_k_changed(self, base_quantity_value):
        self.GenK = base_quantity_value
Here are my doubts:
- I wasn't sure what to specify in

Code: Select all

k = beamsec_obj.GenK.getValueAs("mm^4").Value
as this coefficient is unitless. Should I type something like "mm/mm" or maybe "1" ?
- I_12 must be always zero so I didn't include it in the input parameters, but I don't know if I specified its value properly in the first file
- offset1 and offset2 parameters also have to be zero for this section but I didn't notice them anywhere in FEM module so it's probably not a problem at this point
- two additional changes to the input file have to be made when general section is used (mentioned above: *USER ELEMENT and *ELEMENT, TYPE=U1) - the first one requires adding a whole new keyword after *ELEMENT while the second one requires changing existing parameter from B31 or B32 to U1. Unfortunately, I don't know how to do it yet.
- I hope those are all the necessary files and parts of code that have to be edited to introduce this new functionality
- I'm not sure if the second line of *BEAM SECTION keyword (orientation) is always added in FreeCAD (I hope it is). As I've mentioned, in this case it's required.
- CalculiX also offers box sections, maybe we could add them too.

I will be very grateful for your comments and help with that.
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Arbitrary (general) beam section

Post by bernd »

great start.

How about getting yourself a github account and push these changes to github? It would be very much simpler to pull the changes to my computer and test it.
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Arbitrary (general) beam section

Post by bernd »

have you tested your code? Does it work for you ? I can hardly belive it works ...

Code: Select all

        obj.addProperty(
            "App::PropertyLength",
            "GenIn11",
            "GenBeamSection",
            "set I_11 of the general beam elements"
        )
adds a lengths property ...

Code: Select all

                    in11 = beamsec_obj.GenIn11.getValueAs("mm^4").Value
tries to access the lengths property with a unit mm^4 which is not a lengths. Thus I would expect FreeCAD to return some error. If there is no property for moment of inertia (what I assume) a property fload needs to be added.
User avatar
NewJoker
Veteran
Posts: 3016
Joined: Sun Oct 11, 2020 7:49 pm

Re: Arbitrary (general) beam section

Post by NewJoker »

I haven't tested it yet because I still have some significant doubts that I've listed in my first post. That's why I wanted to discuss it with you first. As you've noticed, the biggest problem is the definition of new variables and constants (I_12 won't be definable for user, it's just 0 at all times). I tried following the syntax of existing beam section definitions but apparently it's not enough when a new variable has a unit of mm^4 instead of regular mm.

I have a github account and I've read the Source code management chapter on wiki so I can try actually implementing it. Your help would be very useful as it's my first approach to FreeCAD code development.
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Arbitrary (general) beam section

Post by bernd »

The reason I have asked for github is because I would like to help ;) This does work best by the use of github. Since you actually do change python only you do not even need to compile.

- get a github account
- fork FreeCAD main repo
- create a new branch
- commit your changes to this branch
- push this branch to your github
- post a link to this branch
User avatar
NewJoker
Veteran
Posts: 3016
Joined: Sun Oct 11, 2020 7:49 pm

Re: Arbitrary (general) beam section

Post by NewJoker »

I added the code from my initial post and pushed the new branch to github. Here's a link for you: https://github.com/FEA-eng/FreeCAD/tree/general_beam

I haven't fixed this issue with unit and property type yet because I wasn't sure where and how to define a new property for moment of inertia. Property for area will probably also have to be changed. And there's that shear coefficient which is unitless, not sure how to handle it.

As I've mentioned before, we will have to add two more features for this to work - new keyword (*USER ELEMENT, ...) and replacement of element type (*ELEMENT, TYPE=U1). There's also one more thing that I forgot to mention before - visualization of those beams in postprocessing. But most likely it won't be a problem since, as far as I know, they are expanded into rectangular section and it should be possible to display results on this representation, maybe even without any additional code edits.

By the way, I've noticed that in the write_femelement_geometry.py file for pipe section the code is:

Code: Select all

section_def = "*BEAM GENERAL SECTION, {}{}{}\n".format(
                        elsetdef,
                        material,
                        section_type
                    )
Shouldn't it be *BEAM SECTION like for other types ? Documentation for CalculiX 2.17 doesn't mention the *BEAM GENERAL SECTION keyword, maybe it's not necessary now and *BEAM SECTION can be used for pipe sections too. Please correct me if I'm wrong.
User avatar
ebrahim raeyat
Posts: 619
Joined: Sun Sep 09, 2018 7:00 pm
Location: Iran
Contact:

Re: Arbitrary (general) beam section

Post by ebrahim raeyat »

NewJoker wrote: Thu Aug 05, 2021 6:09 pm
I haven't fixed this issue with unit and property type yet because I wasn't sure where and how to define a new property for moment of inertia. Property for area will probably also have to be changed. And there's that shear coefficient which is unitless, not sure how to handle it.
As we discussed in pm, I done this part:
1.PNG
1.PNG (13.9 KiB) Viewed 1244 times
2.PNG
2.PNG (20.31 KiB) Viewed 1244 times
NewJoker wrote: Thu Aug 05, 2021 6:09 pm As I've mentioned before, we will have to add two more features for this to work - new keyword (*USER ELEMENT, ...) and replacement of element type (*ELEMENT, TYPE=U1).
again, as we discussed (and brings it here for others to know), writing of ccx is in C++ code in https://github.com/FreeCAD/FreeCAD/blob ... .cpp#L2201.

one dirty solution is that when we have General beam section, in Mod\Fem\femsolver\calculix\write_femelement_geometry.py python code, change the element number type that has same ELSET in General Beam section Block:

Code: Select all

***********************************************************
** Element sets for materials and FEM element type (solid, shell, beam, fluid)
*ELSET,ELSET=M0B0RstdD0
1,

***********************************************************
** Sections
*BEAM SECTION, ELSET=M0B0RstdD0, MATERIAL=MaterialSolid, SECTION=GENERAL
10,11,0,10,10
-3.19744231092E-15, 1, 0[code]
[/code]

@bernd do you have any suggestion? thanks.


example of current development:
FEMMeshGmsh.zip
(1.36 KiB) Downloaded 29 times
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Arbitrary (general) beam section

Post by bernd »

ebrahim raeyat wrote: Thu Nov 03, 2022 1:15 am again, as we discussed (and brings it here for others to know), writing of ccx is in C++ code in https://github.com/FreeCAD/FreeCAD/blob ... .cpp#L2201.

one dirty solution is that when we have General beam section, in Mod\Fem\femsolver\calculix\write_femelement_geometry.py python code, change the element number type that has same ELSET in General Beam section Block:

Code: Select all

***********************************************************
** Element sets for materials and FEM element type (solid, shell, beam, fluid)
*ELSET,ELSET=M0B0RstdD0
1,

***********************************************************
** Sections
*BEAM SECTION, ELSET=M0B0RstdD0, MATERIAL=MaterialSolid, SECTION=GENERAL
10,11,0,10,10
-3.19744231092E-15, 1, 0[code]
[/code]
Only the mesh writing code is in C++ all other ccx solver code is Python. IMHO the problem is not a show stopper. As a fast hack I would read the inp file find the line and replace the element type. I am aware it is a hack, but it would work. If all other stuff of beams are implemented we could fix this. I even might have a look at the C++ if Python beam development would go on.

Here a example of similar problem done in master already ... https://github.com/FreeCAD/FreeCAD/blob ... py#L50-L52 and https://github.com/FreeCAD/FreeCAD/blob ... 2344-L2359
User avatar
ebrahim raeyat
Posts: 619
Joined: Sun Sep 09, 2018 7:00 pm
Location: Iran
Contact:

Re: Arbitrary (general) beam section

Post by ebrahim raeyat »

bernd wrote: Thu Nov 03, 2022 8:53 am Here a example of similar problem done in master already ... https://github.com/FreeCAD/FreeCAD/blob ... py#L50-L52 and https://github.com/FreeCAD/FreeCAD/blob ... 2344-L2359
Thanks. These are very helpful. I will try to fix it.
Post Reply