3 technical questions about how the 6 DOF are implemented

[David_Gallardo]

Hello everyone,

First of all, congratulations to Daniel Piker for the incredible work. I am a structures teacher and I am researching the use of particle systems (PS) as a virtual laboratory for teaching structures.

I have studied the work of Baraff & Witkin (1997), Martini (2001-2011), Kilian (2004), Ochsendorf (2005) and Ramos (2011), which are the fundamentals of the use of particle systems in structures. I understand the convergence problems for very stiff systems (such as building structures), but I believe that FEA can still be used (for example, by making the materials flexible), to visually convey the structural behavior in a qualitative way. I don’t believe that SPs can replace professional FEM programs, but the ability to show the evolution over time of the SPs, together with the direct treatment of non-linear behavior, makes them very interesting for teaching.

Knowing the complexity of the subject, once again, I repeat my congratulations to Daniel Piker for his extraordinary work. The incorporation of the 6 degrees of freedom is one of the greatest contributions to being able to go beyond “form finding” and use it as a structural mechanics laboratory. However, I am trying to understand how it is implemented, and, in the absence of technical documentation, I am having difficulty resolving certain doubts.

I would love to open this thread to gradually incorporate technical information that allows for a more formal and academic use of this excellent tool. To better explain the verification, I have developed a script (Test_kNm.gh) that allows you to obtain the deformation (with a magnification factor to choose from), and the stress laws (axial, shear and moment). Each of the results has a toggle to activate or deactivate its display. The script allows you to change the length of the beam, the number of subdivisions, the load, the dimensions of the cross-section, and the support conditions. I have included explanatory texts where I thought it was necessary to better understand the script.

Test_kNm.gh (108.8 KB)

[David_Gallardo]

So, for @DanielPiker or anyone that could give me some hints, these are the 3 most important questions I have in this example.

• The BEAM component is defined based on the modulus of elasticity of the material E, and the geometric properties of the section (area and moments of inertia). But it is not clear in which units these magnitudes should be entered. I have done tests changing the units of the definition ([kN,m] and [N,mm]), and very different results are obtained (in N and mm, it does not converge), so I understand that they have an influence. To explain it better, the same script in Nmm (test_Nmm.gh). They are identical, except that all the distance and force measurements have different units. In principle, both should converge and give identical results (with the change of units, obviously). But the kNm one converges, and the Nmm one does not. That’s why I suspect there’s something in the solver that depends on the units, and I’d need to know which ones are being used.

Test_Nmm.gh (59.3 KB)

• The results of the SOLVER component can be extracted and analyzed, which is essential to understanding the stresses in the structure. It seems that the bending moments at both ends of each BEAM are obtained (at least that’s what I’ve seen in several videos and tutorials). But the values do not correspond to the moments that can be calculated by hand. I wonder if the values in the results output are actually the curvatures. In structures, bending moments can be obtained from the curvatures by multiplying them by the bending stiffness (E·I). But the values are not correct either. They seem proportional to the moments because the graph is equivalent to the correct moment equation. But the values are not. I don’t know what the values obtained correspond to.
  In fact, if the length of the beam is changed and the calculation is restarted, the same values are obtained every time. In other words, they do not seem to depend on the span of the beam, when in reality the values should depend on the square of the span of the beam.
  However, the values do change when the number of subdivisions of the member is changed. It seems that if the number of subdivisions is halved (for example, from 32 to 16, without changing the span of the beam), the moments are divided by 4.0, and if it is reduced to a quarter (from 32 to 8 divisions), the values obtained are divided by 42.0. The results should vary slightly depending on the size of the discretization, but these huge differences should not occur, because the beam maintains the same span and the same load, so the moments (and the curvatures) should give very similar values.

• The SUPPORT component allows for different types of connections to the outside. I understand that they have been modeled as springs, and the “Strength” parameter indicates the stiffness of the spring (again, the units are in doubt). The default value of 100 is very small compared to the stiffness and loads of real structures. I had to increase it to 1e+12, so that they really are restrained points of movement. However, there must be some problem with the precision limit of floating point numbers, because even though I try to release some rotation, depending on other parameters (such as the load or the rigidity of the BEAM), the rotation is released or not. Is it possible that the release (“False” in the corresponding parameter) does not consider a null rigidity, but a very small one, and that when multiplied by the “Strength” it ends up giving a non-null value? If I change the Strength value to 1e+09, it does work, but with 1e+12, it doesn’t. I guess that if it considers zero rigidity (0.0), it must give a division by zero problem, right? But maybe it could be solved with some conditional for completely freed movements.

I have more questions about technical issues, but I wouldn’t want to overdo it.

I am very grateful for all the work that has been done.

[Tom_P]

some of your questions might be solved by looking at the code of the sample kangaroo 2 Goals on Github
hope this helps - kind regards -tom

[David_Gallardo]

Hi Tom

I already checked it, but i could not find any clues there. In fact, I don’t see there any reference to the BEAM component, or to the SUPPORT component. The most related files are only Anchor.cs or Spring.cs, but i think they correspond to other components, not the 6dof Support.

And the issue with the units and the results of the solver (bending moments) are also not considered in these code files.

Sorry, perhaps I’m not seeing it, and you can give me more details.

Thanks for the help!

[David_Gallardo]

Hello again everyone

Since I wrote the previous post, I have carried out numerous tests and made progress on the issue.

The first thing I want to say is that the work of Cecilie Brandt-Olsen (K2-Engineering) is extraordinary. I have read her thesis and, although there is not much detailed information on the 6-DOF elements, the results obtained are correct, verified with independent software (SAP2000). The data entry is clearly indicated (in N and mm). The output of results, depending on the case, in kN and m, or in N and mm. But it is perfectly indicated.

However, my interest is in understanding why Kangaroo2 (on which K2-Engineering is implemented) does not seem to give correct results. For a qualitative representation, it does work, but, according to the previous example, the values cannot be used from a structural engineering point of view.

Of course, assuming that I have understood correctly how the results are obtained using the BEAM (6DOF) component. As I have seen in different examples, the result provided by the SOLVER for each BEAM component consists of 7 elements:
[0 = initial plane, 1 = final plane, 2 = initial principal moment, 3 = final principal moment, 4 = initial secondary moment, 5 = final secondary moment, 6 = torsional moment]. With the first two, the position of each end in the deformation can be known, so the elastic curve of the deformation can be constructed. The principal bending moments are obtained from the next two elements.

In the simple example in my previous query (a single bar, fixed at both ends, with a uniformly distributed load), the values of elastic deformation and stress laws are easily calculable, considering the material properties of the cross-section, so it is possible to analyze how much the Kangaroo2 result deviates from the correct value.

In fact, the deformation generated by Kangaroo2 and its moment laws do not fit together. The function that represents the elastic deformed curve can be derived twice, and the bending moment law is obtained (well, the curvatures, which multiplied by the bending stiffness EI, are transformed into bending moments). But neither of the two results coincides with the correct solution.

I have found that the deformation obtained in results [0] and [1] is in proportion to (L/N)^3 with respect to the correct solution, where L is the length of the bar and N is the number of discretization divisions. On the other hand, the moment laws obtained in results [2] and [3] are in proportion to (L/N)^2 with respect to the correct solution. However, in both cases there appears to be an additional proportionality constant for which I have found no explanation.

In short, maybe I am extracting the results incorrectly, or entering the data wrongly, and I would like to be able to do it correctly.

In case nobody can explain what is going on, I will use K2-Engineering, but, honestly, I think it would be great to be able to use Kangaroo2 directly with real values of material properties, and get correct results.

Thank you very much!

[arqmorelos1]

Hello @David_Gallardo ,

I have repeatedly read @DanielPiker 's comments recommending to stick to m, kg, N and Pa as units when using kangaroo, so I would personally stick to do the Unit conversion before and after the solver.

Regarding the Nmm not converging there is a Threshold input on the solver which is by default set to 1e-15, I believe this is because normally people would work using meters as units, and it’s much too low a value for your model in mm to converge.

I tried looking at your gh files but you’re wire management makes it very difficult to follow along, but i would think the scripts behaving differently for different units should be about some unit not properly converted.

I did try lowering the threshold on your Nmm file and it does converge, however something is going quite wrong in that file because is very hard for me to visualize anything in Rhino from it.

best regards and good luck
Ricardo

[David_Gallardo]

Thank you @Arqmorelos1.

I agree with you, that the file in Nmm does not make much sense. It is difficult to see, because I’ve been sticking to mm also in the geometry. So the bar is for example 8.000 units in length. And, I also agree that the convergence depends on the threshold, and it is not possible when everything is in mm. So, let’s forget about this model.

With K2-Engineering, the geometry has to be in meters, but all the material and section properties have to be in N and mm (MPa = N/mm2).

I will simplify the gh file (use only N and meter), document it better, and check if then it works correctly. If so, I will inform here asap. If not, I will also give feedback, and post the simplified gh-file.

Thank you a lot for your help! It is really appreciated!

[David_Gallardo]

Hi again, @arqmorelos1

Here is a simplified ghFile.
250303 K2 Test Nm.gh (40.8 KB)

You may use the RemoteControlPanel to play with the different parameters, and see the difference between the expected and the calculated values.

As you may see, the problem is not in the units. I have changed everything to be in N and m, but the results still don’t relate correctly to the correct ones.

It can’t be only a problem in the units, because the solver_results don’t change when changing the length of the beam. The same value is obtained for every beam length, what can’t be correct.

If the length is mantained constant, different results arrive changing the number of divisions. And not in a way that is consequence of the more or less detailed model.

As I mentioned before, I managed to detect that the results of the deflection are proportional to (L/N)^3 with the correct ones, and the results of the moments are proportional to (L/N)^2 with the correct ones. So, if you double the length, and also double the number of divisions, you get the same results from the solver. But if you double the length, being the number of divisions constante, you get 8 times more faktor of proportion between the correct and the calculated deflection, and 4 times more factor of proportion between the correct and the calculated moments.

This does not make any sense. I think only @DanielPiker may have an explanation.

Anyway, thank you for listening!