Построение частотной диаграммы для определения критических скоростей для неуравновешенного ротора в Fidesys 5.1

Интеграция Фидесис с Питоном позволяет проводить обработку результатов, полученных в серии расчетов, с помощью средств визуализации Matplotlib, в частности - построить диаграмму критических скоростей, полученную из серии расчетов собственных частот в модальном анализе с преднагружением угловой скоростью (Prestressed Modal). Необходимо построить модель в Фидесис, а затем перенести ее в скрипт питона и подключить необходимую обработку.

import matplotlib #Библиотека для работы с графиками
matplotlib.use('TKAgg') #Исправляем возможную ошибку по pyQT5
import matplotlib.pyplot as plt #Библиотека для работы с графиками
import matplotlib.ticker as frmt #Библиотека для форматирования графиков
import numpy as np #Библиотека для работы с массивами
import vtk  # Библиотека работы с выходными данными расчетов
from vtk.util.numpy_support import vtk_to_numpy  # Модуль для преобразования результатов
import sys  # Cистемная библиотека
import os  # Cистемная библиотека

fidesys_path = r'C:\Program Files\Fidesys\CAE-Fidesys-5.1'  # Расположение Фидесиса
base_dir = os.path.dirname(os.path.abspath(__file__))  # Директория где лежит скрипт
pvd_file = os.path.join(base_dir, '1.pvd')  # Файл ссылок на результаты
prep_path = os.path.join(fidesys_path, 'preprocessor', 'bin')  # Директория, где препроцессор
os.environ['PATH'] += prep_path  # Добавление пути к препроцессору в PATH
sys.path.append(prep_path)  # Добавление пути к препроцессору в PATH

#Пытаемся импортировать библиотеки Фидесис
try: #Блок попыток
    import cubit                     # Библиотека препроцессинга
    import fidesys                   # Библиотека Фидесиса
except ModuleNotFoundError: #Если не вышло, то выводим в консоль сообщение
    print("В скрипте указан следующий путь к Фидесис: ", fidesys_path)
    print("Укажите в скрипте путь к вашей версии Фидесис. Вероятно ваша версия отличается от указанной или установлена в другую директорию.")
    sys.exit(1)

cubit.init([""])  # Инициализация функций препроцессора
fc = fidesys.FidesysComponent()  # Создание обязательного компонента Фидесис fc
fc.init_application(prep_path) # !!!!!Инициализация для версии 5.1 (для 5.0 и ниже заменить на fc.initApplication(prep_path) )
fc.start_up_no_args()  # Запуск обязательного компонента Фидесис fc

ArrayOfOmega = [1, 30, 50, 70, 75, 85, 100, 150] #Массив частот вращения рад/с
ArrayOfHz = [Om * 0.1592 for Om in ArrayOfOmega] #Массив с пересчитанными частотами в Гц
ArrayOfRPM = [Om * 9.55 for Om in ArrayOfOmega] #Массив с пересчитанными частотами в об/мин
NumOfFreqs = 4 #Число собственных форм для определения
ArrayOfFreqs = [] #Пустой массив для заполнения собственными частотами

for omega in ArrayOfOmega:  #"Пробегаем" в цикле по каждому числу из массива частот вращения
	print("Частота вращения, рад/с: ", omega)  # Пишем в консоль какая частота вращения

	# ---------Начало вставляемого скрипта из Фидесис-------------
	fidesys.cmd('reset')
	fidesys.cmd('create Cylinder height 1 radius 0.01 ')
	fidesys.cmd('create Cylinder height 0.02 radius 0.25 ')
	fidesys.cmd('remove overlap volume 2 1 modify volume 2')
	fidesys.cmd('imprint volume all ')
	fidesys.cmd('merge volume all ')
	fidesys.cmd('volume all size auto factor 3')
	fidesys.cmd('mesh volume all')
	fidesys.cmd('create material 1 from \'Углеродистая сталь\'')
	fidesys.cmd('set duplicate block elements off')
	fidesys.cmd('block 1 add volume all')
	fidesys.cmd('block \'Block 1\' material 1 cs 1 element solid order 1')
	fidesys.cmd('set duplicate block elements off')
	fidesys.cmd('block 2 add node 18102 ')
	fidesys.cmd('create lumpmass properties 1')
	fidesys.cmd('modify lumpmass properties 1 mass 10')
	fidesys.cmd('block \'Block 2\' cs 1 element lumpmass')
	fidesys.cmd('block \'Block 2\' lumpmass properties 1')
	fidesys.cmd('create displacement  on surface 3 2  dof all fix  ')
	fidesys.cmd('create angular velocity global')
	fidesys.cmd('modify angular velocity 1 dof 3 value ' + str(omega))
	fidesys.cmd('analysis type eigenfrequencies findefs elasticity dim3 preload on')

	# ---------Конец вставляемого скрипта из Фидесис-------------

	output_pvd_path = os.path.join(base_dir + "\\" + str(omega) + ".pvd")  # Объявляем директорию и файл сохранения
	print("strarting calculation to " + output_pvd_path)  # Выводим в консоль директорию и файл сохранения
	fidesys.cmd(
		"calculation start path '" + output_pvd_path + "'")  # Просим Фидесис начать расчет в указанную директорию

	print("				  ")
	print("Расчет успешно завершен!")
	
	reader = vtk.vtkXMLUnstructuredGridReader()  # Подключаем читалку
	print("Читаем результаты из ",
		  str(base_dir) + "\\" + str(omega) + r"\case2_step0001_substep0001.vtu")  # Пишет откуда берем результаты
	filename = os.path.join(str(base_dir) + "\\" + str(omega) +r"\case2_step0001_substep0001.vtu")  # Указываем путь к файлу
	reader.SetFileName(filename)  # Подключаем путь к читалке и читаем
	reader.Update()  # Needed because of GetScalarRange
	grid = reader.GetOutput()  # Забираем выходные данные
	point_data = grid.GetPointData()  # Забираем данные для точек

	print("				  ")
	print("Частоты:") # !!!Если в консоль после этой надписи будет выведена какая-то надпись, а не частоты, то это имя массива на данной позиции. Ниже комментарий с указанием где и что менять
	#Делаем функцию, которая "выдергивает" значения частот из результатов
	TempArrayOfFreqs = [] #Создаем временный массив для хранения частот
	for freqNum in range(NumOfFreqs):
		freq = ''
		name = point_data.GetArray(14+freqNum).GetName() # !!!!!Если будет ошибка, то увеличить (или уменьшить) на 1 пока не заработает
		print(name) #Проверяем имя массива, что оно содержит частоту, если нет, то на строке выше меняем цифру
		for ch in (name): #Пробегаем по символам в строке, чтобы вытащить из имени частоту
			if ch.isdigit() or ch == '.': #Если находим число или точку, то ее тоже забираем в число
				freq=freq+str(ch)
			if ch == '(': #Если находим открытую скобку, то опустошаем массив
				freq = ''
		TempArrayOfFreqs.append(float(freq)) #Преобразуем значение СЧ из текста в число и добавляем во временный массив
		print(freq)
	ArrayOfFreqs.append(TempArrayOfFreqs) #Из временного массива берем частоту и добавляем в основной массив
	print("				  ")
	#Здесь цикл закрывается
	
fc.delete_application() # !!!Команда удаления задачи из памяти для версии 5.1 (для 5.0 и ниже заменить на fc.deleteApplication())
print("Готово!")

#Создаем объект графика и осей
fig = plt.figure()
axs = fig.add_subplot()

ArrayOfFreqsT = np.transpose(ArrayOfFreqs) #Меняем строки со столбцами в массиве с наборами частот
axs.xaxis.set_major_locator(frmt.FixedLocator(ArrayOfHz)) #Указываем на каких позициях по частотам в Гц будут подписи в нужных единицах
axs.xaxis.set_major_formatter(frmt.FixedFormatter(ArrayOfOmega)) #Задаем подписи для вышеуказанных позиций в рад/с
#axs.xaxis.set_major_formatter(frmt.FixedFormatter(ArrayOfRPM)) #Задаем подписи для вышеуказанных позиций в об/мин

#Рисуем графики изменения собственных частот от частот вращения
freqN = 1
for freq in ArrayOfFreqsT:
	plt.plot(ArrayOfHz,freq,label=str(freqN)+' частота')
	freqN+=1
#Рисуем линию совпадения частот вращения, пересчитанных в Гц с собственными
plt.plot(ArrayOfHz,ArrayOfHz,label='Критическая линия',color = 'r')
#plt.plot(ArrayOfHz,[Hz * 2 for Hz in ArrayOfHz],label='Критическая линия 2x')
	
#Делаем подписи осей для графика
plt.title("Диаграмма критических скоростей")
plt.xlabel("Частота вращения, рад/с")
#plt.xlabel("Частота вращения, об/мин")
plt.ylabel("Собственная частота, Гц")

#Рисуем сетку, легенду и сам график
plt.grid()
plt.legend()
plt.show()

Результат будет иметь следующий вид: