byuflowlab / GXBeam.jl

Pure Julia Implementation of Geometrically Exact Beam Theory
MIT License
87 stars 18 forks source link

Rotation angle #69

Closed luizpancini closed 2 years ago

luizpancini commented 2 years ago

Hi,

First I'd like to say that I really appreciate the work you guys put into this package, great stuff, great initiative on making it open source.

While comparing the results with some of the literature cases on geometrically exact beams, I found that the nodal rotation angles were a little off. For instance, in the example Nonlinear Analysis of a Cantilever Subjected to a Constant Moment, the rotation angle theta_y should vary linearly with the load factor lambda, but it doesn't in the program.

After some research, realized that the rotation parameters, and not the rotation angles, were being output in the functions extract_point_states! and extract_element_state., check e.g. equation 10 of your Ref. [2].

The solution is simple though, just make phi = 4*atan(theta/4), and output phi instead of theta as the rotation angles of the cross-section. In the attached figure, phi is the rotation angle and theta is the rotation parameter for the aforementioned example.

I tried to make the correction on a fork and set a pull request, but I am apparently too much of a noob in Julia even to do just that!

[2] Wang, Q., & Yu, W. (2017). Geometrically nonlinear analysis of composite beams using Wiener-Milenković parameters. Journal of Renewable and Sustainable Energy, 9(3), 033306.

example

taylormcd commented 2 years ago

I'm glad you like the package. Converting from rotation parameters to angles makes sense for the 2D case because the axis of rotation is obvious, but the same isn't true for the general case. Converting to Euler angles in the 3D case would require choosing a rotation convention. I feel like its more natural to just construct a rotation matrix from the Wiener-Milenkovic parameters directly rather than convert to some other rotation parameterization first. Is there somewhere in the documentation where I said angles when I should have said rotation parameters?

luizpancini commented 2 years ago

The documentation is clear about the outputs being the rotation parameters, it was just my bias that made me think they were the angles. However, there is no conversion of parametrization involved in making that step, it directly gives the rotation angles in the global, body-attached, "a" frame.

taylormcd commented 2 years ago

The rotation parameters are defined by c_i = 4 tan(phi/4) n_i where c is the Wiener Milenkovic rotation parameters, phi is the rotation angle, and n is the rotation axis. The rotation axis n needs to be known in order to directly solve for the angle using this equation.

taylormcd commented 2 years ago

It isn't exactly the most intuitive of rotation conventions so perhaps the user interface should instead accept any of the rotation parameterizations provided by the Rotations.jl package and convert to Wiener Milenkovic parameters internally.

luizpancini commented 2 years ago

Yeah, providing that option would make it easier to define a complicated assembly with several beams going in several directions, and easier to interpret the results. What I meant before was that you may get the parameters for any geometrical parametrization (e.g. Euler angles) directly from the rotation tensor.