ansys / pymapdl-reader

Legacy binary interface to MAPDL binary files.
https://reader.docs.pyansys.com
MIT License
47 stars 23 forks source link

Cyclic analysis: forward/backward mode direction is inconsistent #86

Open mjc865 opened 2 years ago

mjc865 commented 2 years ago

For a cyclic analysis, it looks like there's some inconsistency between the harmonic index of the result set, and the harmonic index as calculated from the sector boundary. Sometimes they match, and sometimes they don't.

Here's an uploaded results file (file_cyclic_test.rst) generated from Mechanical 2020R2: http://filedropper.com/325u3Hvx (simple model with 36 base sector nodes, 8 sectors, with results calculated at 1 ND)

And some code to calculate the harmonic index from two nodes on the sector boundary, and compare with the value from the result set:

from ansys.mapdl import reader as pymapdl_reader
import numpy as np
from math import pi

rstfile = 'D:/temp/file_cyclic_test.rst'  # replace with path to uploaded RST file

result = pymapdl_reader.read_binary(rstfile)

low_node_index = 28  # index of low side node
high_node_index = 35  # index of matching high side node

for i_result in range(0, result.n_results):
    node_deflection = result.nodal_displacement(i_result, as_complex=True)

    # only need one component to calculate phase angle
    deflection_low = node_deflection[1][low_node_index, 2]
    deflection_high = node_deflection[1][high_node_index, 2]
    angle_low = np.angle(deflection_low)
    angle_high = np.angle(deflection_high)

    # calculated nodal diameter from sector boundary phase angles
    nd_calculated = (angle_high - angle_low)*result.n_sector / (2*pi)

    # make sure this lies between -n_sector/2 and n_sector/2
    if nd_calculated > result.n_sector/2:
        nd_calculated -= result.n_sector
    if nd_calculated <= -1*result.n_sector/2:
        nd_calculated += result.n_sector

    # round to get rid of small floating point errors
    nd_calculated = round(nd_calculated, 3)

    print('\nResult set ' + str(i_result))
    print('Harmonic index from result set: ' + str(result.harmonic_indices[i_result]))
    print('Harmonic index calculated from boundary: ' + str(nd_calculated))
    if result.harmonic_indices[i_result] == int(nd_calculated):
        print('Match!')
    else:
        print('No match!')

This gives the following result:

Result set 0 Harmonic index from result set: 1 Harmonic index calculated from boundary: 1.0 Match!

Result set 1 Harmonic index from result set: -1 Harmonic index calculated from boundary: -1.0 Match!

Result set 2 Harmonic index from result set: 1 Harmonic index calculated from boundary: 1.0 Match!

Result set 3 Harmonic index from result set: -1 Harmonic index calculated from boundary: -1.0 Match!

Result set 4 Harmonic index from result set: 1 Harmonic index calculated from boundary: -1.0 No match!

Result set 5 Harmonic index from result set: -1 Harmonic index calculated from boundary: 1.0 No match!

Result set 6 Harmonic index from result set: 1 Harmonic index calculated from boundary: 1.0 Match!

Result set 7 Harmonic index from result set: -1 Harmonic index calculated from boundary: -1.0 Match!

Result set 8 Harmonic index from result set: 1 Harmonic index calculated from boundary: -1.0 No match!

Result set 9 Harmonic index from result set: -1 Harmonic index calculated from boundary: 1.0 No match!

Result set 10 Harmonic index from result set: 1 Harmonic index calculated from boundary: 1.0 Match!

Result set 11 Harmonic index from result set: -1 Harmonic index calculated from boundary: -1.0 Match!

Process finished with exit code 0

In the first two mode pairs (0,1 and 2,3), the direction matches between the result set and the boundary calculation, but for the third mode pair (4,5) the direction doesn't.

akaszynski commented 2 years ago

@mjc865, I was away from this repo over can you please upload again (ideally here in the issue)?

mjc865 commented 2 years ago

file_cyclic_test.zip

I was actually looking into this yesterday - it looks like the signs of the harmonic indices aren't stored in the actual RST file, just interpreted from the code:

(rst.py line 220)

# ansys 15 doesn't track negative harmonic indices
    if not np.any(hindex < -1):
        # check if duplicate frequencies
        tvalues = resultheader['time_values']
        for i in range(tvalues.size - 1):
            # adjust tolarance(?)
            if np.isclose(tvalues[i], tvalues[i + 1]):  
                hindex[i + 1] *= -1

I guess you'd have to find matching low/high side nodes and calculate the direction from that, since the direction in ANSYS is inconsistent (first mode of each pair is not always in the positive direction). Workbench automatically creates components for the low/high side nodes but I'm not sure if this is true if the analysis is not run through Workbench.

Also, you may want to check the harmonic index of the result set and exclude the N/2 modes, as these frequencies do not occur in pairs so even if they're close in frequency, they are distinct modes.

akaszynski commented 2 years ago

it looks like the signs of the harmonic indices aren't stored in the actual RST file, just interpreted from the code

That's my understanding as well, and it forced me to "guess" within rst.py.

Components are now written to rst files, so we can check if they exist as MAPDL will store these components when running CYCLIC.