imcs-compsim / meshpy

A general purpose 3D beam finite element input generator
Other
9 stars 3 forks source link

How to calculate the twist angle in create_beam_mesh_helix #87

Closed knarfnitram closed 1 week ago

knarfnitram commented 3 weeks ago

I am currently trying to create a flow diverter with create_beam_mesh_helix. I am therefore creating two helixes in different directions (see variable clockwise), but I am having troubles/ problems understanding the parameter twist_angle. From my understanding, if I choose the twist angle with respect to 2*pi and the length of the helix, I should end up at the same starting point in x and y-direction, but with an offset according to the length of the helix.

However, if we take look at the picture, there's a gap between the two different endpoints (which should coincide from my understanding) Screenshot from 2024-07-12 17-33-41

Code Example: L=1 for clockwise in [-1.0,1.0]: beam_set = create_beam_mesh_helix( mesh_flow_diverter, beam_object, mat, axis=[0.0, 0.0, clockwise], axis_point=[0.0, 0.0, 0.0], start_point=[1.0, 0.0, 0.0], twist_angle=-clockwise*L/(2*np.pi), height_helix=L, n_el=100, )

Am i using it wrong or is there an inaccuray?

davidrudlstorfer commented 3 weeks ago

Thanks for raising the issue - I'm going to look into it ASAP and will add a comment once I found out more

davidrudlstorfer commented 3 weeks ago

Regarding your issue with the helix computation I've found out where you misinterpreted the twist angle:

The angle is the actual angle of the helix, so you need to calculate the angle with the arc tangent as shown in the sketch below:

IMG_2461

For easier reference I've renamed the height of the helix to H and introduced the radius R = 1

With the following code:

import numpy as np

from meshpy import (
    InputFile,
    MaterialReissner,
    Beam3rHerm2Line3,
)

from meshpy.mesh_creation_functions.beam_basic_geometry import create_beam_mesh_helix

def main():

    # define input file
    input_file = InputFile()

    # define material
    mat = MaterialReissner(youngs_modulus=1000, radius=0.02, shear_correction=1.0)

    # height of the helix
    H = 1

    # radius of the helix
    R = 1.0

    for clockwise in [-1.0, 1.0]:
        create_beam_mesh_helix(
            input_file,
            Beam3rHerm2Line3,
            mat,
            axis=[0.0, 0.0, clockwise],
            axis_point=[0.0, 0.0, 0.0],
            start_point=[R, 0.0, 0.0],
            twist_angle=-clockwise * np.arctan(H / (2 * np.pi * R)),
            height_helix=H,
            n_el=100,
        )

    input_file.display_pyvista(
        beam_nodes=True,
        beam_tube=True,
        beam_cross_section_directors=True,
        resolution=1,
        parallel_projection=True,
    )

if __name__ == "__main__":

    main()

I get the correct (and probably what you want) helix as follows:

Screenshot from 2024-07-15 09-51-05

I hope this helps with your geometry creation!

EDIT: with H = 1 and R = 1, $\frac{1}{2 \pi}$ and $arctan(\frac{1}{2 \pi})$ are really close, that's why only a small gap was present with your geometry

knarfnitram commented 3 weeks ago

@davidrudlstorfer, thanks for fast response and the explanation!!

isteinbrecher commented 2 weeks ago

Thanks @davidrudlstorfer for the explanation. Would it make sense to rename twist_angle to pitch_angle to avoid confusion?

davidrudlstorfer commented 2 weeks ago

@isteinbrecher I would be open for a renaming if it improves the understanding. In my opinion it is hard to come up with a name that is self explanatory. During the initial research of the helix quantities a few names for the angle came up including:

All of them are used interchangeable (some very random examples: helix angle, twist angle, pitch angle) The most prominent name is "helix angle". In my opinion pitch angle could be confusing as pitch is something different as an angle.

We could rename it to helix angle and also add the synonyms, i.e., pitch angle & twist angle? If you prefer pitch angle I would also be more than happy to rename it to that.

isteinbrecher commented 2 weeks ago

Look like we had a look at the same examples. Another idea, we could leave the names as it is, but make twist_angle also optional and require two of the three parameters (twist_angle, helix_height, turns) to be set. Then we can do the mentioned conversion in this thread directly inside the create_helix function. What do you think?

davidrudlstorfer commented 1 week ago

@knarfnitram maybe it helps with your setup

With #93, the above discussed mesh creation simplifies to:

import numpy as np

from meshpy import (
    InputFile,
    MaterialReissner,
    Beam3rHerm2Line3,
)

from meshpy.mesh_creation_functions.beam_basic_geometry import create_beam_mesh_helix

def main():

    # define input file
    input_file = InputFile()

    # define material
    mat = MaterialReissner(youngs_modulus=1000, radius=0.02, shear_correction=1.0)

    for turn in [-1.0, 1.0]:
        create_beam_mesh_helix(
            input_file,
            Beam3rHerm2Line3,
            mat,
            axis=[0.0, 0.0, 1.0],
            axis_point=[0.0, 0.0, 0.0],
            start_point=[1.0, 0.0, 0.0],
            height_helix=1.0,
            turns=turn,
            n_el=100,
        )

    input_file.display_pyvista(
        beam_nodes=True,
        beam_tube=True,
        beam_cross_section_directors=True,
        resolution=1,
        parallel_projection=True,
    )

if __name__ == "__main__":

    main()
knarfnitram commented 1 week ago

hey @davidrudlstorfer, thanks for adjusting the creation of a helix with beams, I think this clarifies and simplifies a lot!