fidesys
В данной статье показан пример того, как можно создать пользовательскую утилиту с графическим интерфейсом для расчета критических скоростей вращения системы вал-ротор с учетом гироскопического эффекта и построения диаграммы Кэмпбелла.
Для корректной работы python-файла, сначало нужно установить библиотеки pyQt5=5.15.1 (именно 5.15.1, поскольку выше версии предпроцессор фидесиса не поддерживает), numpy, scipy, meshio, pyqtgraph: pip install pyQt5==5.15.1, pip install numpy, pip install scipy, pip install meshio, pip install pyqtgraph.
Создайте новый файл python. Для этого запустите среду разработки IDLE из меню «Пуск» или папки, где установлен Python. Выберите File - New File.
Cкопируйте в открывшееся пустое окно следующий скрипт:
import shutil,os
JOIN = os.path.join
# function clear working folder
def clear_working_folder(dir,file_name):
path_working_folder = os.listdir(dir)
for element in path_working_folder:
element_path = JOIN(dir, element)
if element.startswith(file_name):
if os.path.isfile(element_path):
os.remove(element_path)
if os.path.isdir(element_path):
shutil.rmtree(element_path)
else:
continue
Сохраните его с названием clear_working_folder.py
Затем создайте новый файл python в той же папке, где вы сохранили предыдущий файл:
# VM 247
import numpy as np
import os,sys,time
from clear_working_folder import JOIN
FIDESYS_PREP_DIR = 'C:/Program Files/Fidesys/CAE-Fidesys-6.0/preprocessor/bin'
prep_paths = [
FIDESYS_PREP_DIR,
JOIN(FIDESYS_PREP_DIR, 'plugins'),
JOIN(FIDESYS_PREP_DIR, 'acis', 'code', 'bin')
] # Preprocessor directory
calc_path = JOIN(FIDESYS_PREP_DIR, 'FidesysCalc.exe') # solution directory
for path in prep_paths:
os.environ['PATH'] += os.pathsep + path # ADD directories to PATH.
os.add_dll_directory(path)
sys.path.append(path) # ADD directories to PATH.
try:
import cubit, fidesys
except ModuleNotFoundError:
print(f'Модули fydesys, cubit не найдены в директории {FIDESYS_PREP_DIR}')
sys.exit(1)
BASE_DIR = os.getcwd()
RESULT_FILENAME = 'modal'
modes_count = 6
# diametres of circle sections
section_diametr_array = [
1.02E-2,
2.04E-2,
1.52E-2,
4.06E-2,
4.06E-2,
6.6E-2,
6.6E-2,
5.08E-2,
5.08E-2,
2.54E-2,
2.54E-2,
3.04E-2,
3.04E-2,
2.54E-2,
2.54E-2,
7.62E-2,
4.06E-2,
4.06E-2
]
# block ids with pipe section
tube_section_ids = [7,8,18]
keypoints_beam_coords = [
0,
1.27E-2,
5.08E-2,
7.62E-2,
8.89E-2,
10.16E-2,
10.67E-2,
11.43E-2,
12.7E-2,
13.46E-2,
16.51E-2,
19.05E-2,
22.86E-2,
26.67E-2,
28.7E-2,
30.48E-2,
31.5E-2,
34.54E-2,
35.5E-2
]
keypoints_spring_coords = [
16.51E-2,16.51E-2,28.7E-2,28.7E-2
]
def create_beam_blocks():
fidesys.cmd('create material 1')
fidesys.cmd('modify material 1 name \'stal\'')
fidesys.cmd('modify material 1 set property \'MODULUS\' value 2.078e+11')
fidesys.cmd('modify material 1 set property \'SHEAR_MODULUS\' value 1e+12')
fidesys.cmd('modify material 1 set property \'DENSITY\' value 7806')
for keypoint_coord in keypoints_beam_coords:
fidesys.cmd(f'create vertex {keypoint_coord} 0 0')
for keypoint_id in range(1,len(keypoints_beam_coords)+1):
if keypoint_id < len(keypoints_beam_coords):
fidesys.cmd(f'create curve vertex {keypoint_id} {keypoint_id+1}')
else:
break
fidesys.cmd('merge all')
fidesys.cmd('curve all interval 1')
fidesys.cmd('curve all scheme equal')
fidesys.cmd('mesh curve all')
fidesys.cmd('merge all')
for index in range(len(section_diametr_array)):
diametr = section_diametr_array[index]
fidesys.cmd(f'block {index + 1} add curve {index + 1}')
fidesys.cmd(f'create beam properties {index + 1}')
fidesys.cmd(f'modify beam properties {index + 1} type \'Ellipse\'')
fidesys.cmd(f'modify beam properties {index + 1} angle 0.0')
fidesys.cmd(f'modify beam properties {index + 1} ey 0.0')
fidesys.cmd(f'modify beam properties {index + 1} ez 0.0')
fidesys.cmd(f'modify beam properties {index + 1} geom_a {diametr}')
fidesys.cmd(f'modify beam properties {index + 1} geom_b {diametr}')
fidesys.cmd(f'modify beam properties {index + 1} mesh_quality 6')
fidesys.cmd(f'modify beam properties {index + 1} warping_dof off')
fidesys.cmd(f'block {index + 1} material 1 cs 1 element beam order 1')
fidesys.cmd(f'block {index + 1} beam properties {index + 1}')
# modify pipe sections
for tube_section_id in tube_section_ids:
internal_diametr = 0.0304
fidesys.cmd(f'modify beam properties {tube_section_id} type \'Circle With Offset Hole\'')
fidesys.cmd(f'modify beam properties {tube_section_id} angle 0')
fidesys.cmd(f'modify beam properties {tube_section_id} ey 0')
fidesys.cmd(f'modify beam properties {tube_section_id} ez 0')
fidesys.cmd(f'modify beam properties {tube_section_id} geom_d1 {section_diametr_array[tube_section_id - 1]}')
internal_diametr = 0.0356 if tube_section_id == 8 else internal_diametr
fidesys.cmd(f'modify beam properties {tube_section_id} geom_d2 {internal_diametr}')
fidesys.cmd(f'modify beam properties {tube_section_id} geom_e 0')
fidesys.cmd(f'modify beam properties {tube_section_id} mesh_quality 6')
fidesys.cmd(f'modify beam properties {tube_section_id} warping_dof off')
fidesys.cmd('create displacement on vertex all dof 1 dof 4 fix 0')
def create_combin_blocks():
increment = 0
length = len(keypoints_beam_coords)
for keypoints_spring_coord_index in range(len(keypoints_spring_coords)):
if keypoints_spring_coord_index % 2 == 0:
fidesys.cmd(f'create curve location {keypoints_spring_coords[keypoints_spring_coord_index]} 0 0 location {keypoints_spring_coords[keypoints_spring_coord_index]} 1e-03 0')
else:
fidesys.cmd(f'create curve location {keypoints_spring_coords[keypoints_spring_coord_index]} 0 0 location {keypoints_spring_coords[keypoints_spring_coord_index]} 0 1e-03')
increment += 1
combin_line_id = length - 1 + increment
fidesys.cmd('curve all interval 1')
fidesys.cmd('curve all scheme equal')
fidesys.cmd(f'mesh curve {combin_line_id}')
fidesys.cmd(f'block {length } add curve {combin_line_id}')
fidesys.cmd('merge all')
fidesys.cmd('create spring properties 1')
fidesys.cmd('modify spring properties 1 type \'linear_spring\'')
fidesys.cmd('modify spring properties 1 stiffness 4.37e7')
fidesys.cmd('modify spring properties 1 spring_constant_damping 0')
fidesys.cmd('modify spring properties 1 spring_linear_damping 0')
fidesys.cmd('modify spring properties 1 spring_mass 0')
fidesys.cmd(f'block {length} element spring')
fidesys.cmd(f'block {length} spring properties 1')
fidesys.cmd('create displacement on vertex 42 38 44 40 dof all fix')
def create_point_mass_block():
length = len(keypoints_beam_coords)
fidesys.cmd(f'block { length + 1 } add vertex 5')
fidesys.cmd('create lumpmass properties 1')
fidesys.cmd('modify lumpmass properties 1 mass_x 1.401')
fidesys.cmd('modify lumpmass properties 1 mass_y 1.401')
fidesys.cmd('modify lumpmass properties 1 mass_z 1.401')
fidesys.cmd('modify lumpmass properties 1 mass_inertia_x 0.002')
fidesys.cmd('modify lumpmass properties 1 mass_inertia_y 0.00136')
fidesys.cmd('modify lumpmass properties 1 mass_inertia_z 0.00136')
fidesys.cmd(f'block { length + 1 } cs 1 element lumpmass')
fidesys.cmd(f'block { length + 1 } lumpmass properties 1')
def calculate(omega):
modal_result_path = f'{BASE_DIR}/{RESULT_FILENAME}_{round(omega/0.1047,3) }'
cubit.init([''])
fc = fidesys.FidesysComponent()
fc.init_application(FIDESYS_PREP_DIR) # initialization
fc.start_up_no_args() # run fidesys component fc
fidesys.cmd('reset')
create_beam_blocks()
create_combin_blocks()
create_point_mass_block()
fidesys.cmd('create angular velocity global')
fidesys.cmd(f'modify angular velocity 1 dof 1 value {omega}')
fidesys.cmd('analysis type eigenfrequencies dim3 preload off')
fidesys.cmd(f'eigenvalue find {modes_count} smallest')
fidesys.cmd('damping off mass_matrix 0 stiffness_matrix 0')
fidesys.cmd('coriolis on')
fidesys.cmd(f"calculation start path '{modal_result_path}.pvd'")
fc.delete_application()
Сохраните файл с названием calculate_modal.py
Еще раз создайте новый файл python в той же папке, где вы сохранили предыдущий файл:
from PyQt5 import QtCore,QtTest
from PyQt5.QtCore import QRegExp
from PyQt5.QtGui import QRegExpValidator
from PyQt5.QtWidgets import QApplication, QMainWindow,QPushButton,QTableWidget,QTableWidgetItem,QWidget,QGridLayout,QLineEdit,QLabel,QMessageBox,QPlainTextEdit
import sys,os,meshio
import numpy as np
import pyqtgraph as pg
from clear_working_folder import clear_working_folder
from calculate_modal import modes_count,BASE_DIR,RESULT_FILENAME,calculate
from scipy.interpolate import interp1d
from scipy.linalg import solve
from scipy.optimize import root
from random import randrange
class Window(QMainWindow):
rotational_velocities_rpm = np.array([0,35000,70000,105000])
rotational_velocities_radian = rotational_velocities_rpm * 0.1047
rotational_velocities_hz = rotational_velocities_rpm * 0.017
def __init__(self):
super().__init__()
self.setupUi() # Create Widgets
def setupUi(self):
self.setWindowTitle('Campell diagram')
self.setFixedHeight(1000)
self.setFixedWidth(1000)
self.pre_processor_title_create()
self.btn_calculate_create()
self.table_list_rotational_velocities_create()
self.post_processor_title_create()
self.graph_Widget_create()
self.table_list_critical_speeds_create()
self.text_box_ratio_edit_create()
self.text_title_ratio_create()
self.button_replot_graph_and_recalculate_critical_speeds_create()
self.terminal_create()
self.toggle_visibility_widgets()
# empty Function
def empty_function(self,omega):
pass
# BLOCK OF CREATING WIDGETS BEGINS!!
# create table list of rotational velocities
def table_list_rotational_velocities_create(self):
length = len(self.rotational_velocities_rpm)
self.table_list_rotational_velocities = QTableWidget(self)
self.table_list_rotational_velocities.setGeometry(50,50,220,120)
self.table_list_rotational_velocities.setRowCount(length)
self.table_list_rotational_velocities.setColumnCount(1)
self.table_list_rotational_velocities.setHorizontalHeaderLabels(['Rotational velocity (rpm)'])
self.table_list_rotational_velocities.horizontalHeader().setDefaultSectionSize(200)
self.table_list_rotational_velocities.verticalHeader().setDefaultSectionSize(15)
self.insert_rotational_velocities_values_into_table()
#create preprocessor title
def pre_processor_title_create(self):
self.pre_processor_title = QLabel(self)
self.pre_processor_title.setText('PREPROCESSOR (4 VELOCITIES -> 4 STEPS)')
self.pre_processor_title.setStyleSheet('font-size:30px')
self.pre_processor_title.setGeometry(200, 0, 700, 40)
# create button calculate
def btn_calculate_create(self):
self.btn_calculate = QPushButton('Calculate fidesys script',self)
self.btn_calculate.move(50,350)
self.btn_calculate.setGeometry(500, 80, 200, 40)
self.btn_calculate.clicked.connect(self.btn_calculate_handler)
# create postprocessor title
def post_processor_title_create(self):
self.post_processor_title = QLabel(self)
self.post_processor_title.setText('POSTPROCESSOR')
self.post_processor_title.setStyleSheet('font-size:30px')
self.post_processor_title.setGeometry(400, 380, 300, 40)
# create graph widget
def graph_Widget_create(self):
self.plot_graph = pg.PlotWidget()
self.plot_graph.setBackground('w')
self.plot_graph.addLegend(offset=(-650, 5))
self.plot_graph.setStyleSheet("border : none;")
self.widget = QWidget(self)
self.layout = QGridLayout()
self.layout.addWidget(self.plot_graph)
self.widget.setLayout(self.layout)
self.widget.setStyleSheet("border : 5px solid blue;")
self.widget.setGeometry(100,420, 800, 400)
# create table list of critical speeds
def table_list_critical_speeds_create(self):
self.table_list_critical_speeds = QTableWidget(self)
self.table_list_critical_speeds.setGeometry(50,850,220,120)
self.table_list_critical_speeds.setColumnCount(1)
self.table_list_critical_speeds.setHorizontalHeaderLabels(['Critical speeds (rpm)'])
self.table_list_critical_speeds.horizontalHeader().setDefaultSectionSize(200)
self.table_list_critical_speeds.verticalHeader().setDefaultSectionSize(15)
# create text box edit ratio of critical line
def text_box_ratio_edit_create(self):
self.text_box_ratio_edit = QLineEdit(self)
self.text_box_ratio_edit.setGeometry(450,925,50,20)
self.text_box_ratio_edit.setAlignment(QtCore.Qt.AlignCenter)
validator = QRegExpValidator(QRegExp("\d+\.?\d*e?-?\d+"))
self.text_box_ratio_edit.setValidator(validator)
self.text_box_ratio_edit.setText('1')
# create button replot campell's diagram and recalculate critical speeds
def button_replot_graph_and_recalculate_critical_speeds_create(self):
self.button_replot_graph_and_recalculate_critical_speeds = QPushButton('Replot',self)
self.button_replot_graph_and_recalculate_critical_speeds.setGeometry(550,920,150,30)
self.button_replot_graph_and_recalculate_critical_speeds.clicked.connect(self.button_replot_graph_and_recalculate_critical_speeds_handler)
# create text title 'RATIO'
def text_title_ratio_create(self):
self.text_title_ratio = QLabel(self)
self.text_title_ratio.setGeometry(400,925,50,20)
self.text_title_ratio.setText('RATIO')
# create terminal
def terminal_create(self):
self.terminal = QPlainTextEdit(self)
self.terminal.setGeometry(100,180, 800, 200)
# BLOCK OF CREATING WIDGETS FINISH!!
# control visibility of Widgets
def toggle_visibility_widgets(self,is_visible = False):
self.button_replot_graph_and_recalculate_critical_speeds.setVisible(is_visible)
self.text_box_ratio_edit.setVisible(is_visible)
self.text_title_ratio.setVisible(is_visible)
self.table_list_critical_speeds.setVisible(is_visible)
# plot diagram campell
def plot_campell_diagramm(self,modes_dict,ratio = 1):
def create_random_color(): # random hex color
symbols = [1,2,3,4,5,6,7,8,9,'A','B','C','D','E']
color = '#'
for i in range(0,6):
random_index = randrange(0,len(symbols))
random_symbol = symbols[random_index]
color = color + str(random_symbol)
return color
self.plot_graph.clear()
for mode_index in range(1,modes_count+1):
line_color = create_random_color()
self.plot_graph.plot(self.rotational_velocities_rpm,modes_dict[f'{mode_index} mode'], fillLevel = 10,pen = pg.mkPen(color = line_color,width = 3),name = f'Mode_{mode_index}')
self.plot_graph.plot(self.rotational_velocities_rpm,self.rotational_velocities_hz/float(ratio),pen = pg.mkPen(color = "black",width = 4))
# insert_rotational_velocities_values_into_table
def insert_rotational_velocities_values_into_table(self):
table = self.table_list_rotational_velocities
length = len(self.rotational_velocities_rpm)
for index in range(length):
table.setItem(index,0,QTableWidgetItem(str(self.rotational_velocities_rpm[index])))
table.item(index,0).setFlags(QtCore.Qt.ItemIsEditable)
table.item(index,0).setTextAlignment(QtCore.Qt.AlignHCenter)
# get damping frequencies from vtu
def get_frequencies_from_vtu(self,calculate_fidesys_script_function):
modes_dict = {}
modes_dict_for_plotting = {}
for omega in self.rotational_velocities_radian:
modal_result_path_result_vtu = f'{BASE_DIR}/{RESULT_FILENAME}_{round(omega/0.1047,3)}/case1_step0001_substep0001.vtu'
calculate_fidesys_script_function(omega)
try:
mesh = meshio.read(f'{modal_result_path_result_vtu}')
except meshio._exceptions.ReadError:
self.terminal.insertPlainText(f'Error!!!.Vtu file from {modal_result_path_result_vtu} is not Found\n')
clear_working_folder(BASE_DIR,RESULT_FILENAME)
self.btn_calculate.setEnabled(True)
raise meshio._exceptions.ReadError('Ошибка! Vtu-Файл не найден!')
damping_freqs = mesh.field_data['Eigen Values Extended']
damping_freqs = list(damping_freqs[0])
modes_dict[f'{omega} rad/s'] = damping_freqs
for mode_index in range(1,modes_count+1):
modes_dict_for_plotting[f'{mode_index} mode'] = []
for key in modes_dict:
element = round(modes_dict[key][mode_index-1],3)
modes_dict_for_plotting[f'{mode_index} mode'].append(element)
for key,value in modes_dict_for_plotting.items():
modes_dict_for_plotting[key] = np.array(value)
return modes_dict_for_plotting
# calculate damping Frequencies
def calculate_damping_frequencies(self,omega):
self.terminal.insertPlainText(f'Calculation with {round(omega/0.1047,3)} rpm begin\n')
calculate(omega)
QtTest.QTest.qWait(3000)
self.terminal.insertPlainText(f'Calculation with {round(omega/0.1047,3)} rpm finished succesfully\n')
# calculate critical speeds
def critital_speeds_calculate(self,modes_dict,ratio = 1):
self.toggle_visibility_widgets(True)
critical_line_matrix = np.array([[self.rotational_velocities_rpm[0],1],[self.rotational_velocities_rpm[-1],1]])
critical_line_vector = np.array([self.rotational_velocities_hz[0]/float(ratio),self.rotational_velocities_hz[-1]/float(ratio)])
[k1,b1] = solve(critical_line_matrix,critical_line_vector) # equation of line y = k1*x+ b1. Find coeffs k1,b1
def intersect_function(x): # point of intersection of two linear functions
return k*x + b - k1*x - b1
critical_speeds = []
for key in modes_dict:
array = modes_dict[key]
for index in range(0,len(array)-1):
current_rotational_velocity = self.rotational_velocities_rpm[index]
next_rotational_velocity = self.rotational_velocities_rpm[index + 1]
current_speed = array[index]
next_speed = array[index + 1]
interval_rotational_velocity = np.array([current_rotational_velocity,next_rotational_velocity])
interval_speed = np.array([current_speed,next_speed])
linear_interpolation =interp1d(interval_rotational_velocity,interval_speed,kind = 'linear')
matrix = np.array([[current_rotational_velocity,1],[next_rotational_velocity,1]])
vector = np.array([current_speed,next_speed])
[k,b] = solve(matrix,vector) # equation of line y = k*x+b. Find coeffs k,b
sol = root(intersect_function, x0=0) # find intersect of two functions in current_interval
if current_rotational_velocity <= sol.x[0] <= next_rotational_velocity:
critical_speed = round(float(linear_interpolation(sol.x[0])),3)
critical_speeds.append(critical_speed)
else:
continue
critical_speeds = np.array(critical_speeds)/0.017 # critical speeds in rpm
return np.round(critical_speeds,3)
# clear critical speeds table
def table_list_critical_speeds_clear(self):
table = self.table_list_critical_speeds
while (table.rowCount() > 0):
table.removeRow(0)
# insert critical speeds into table
def insert_critical_speeds_into_table(self,modes_dict,ratio = 1):
self.table_list_critical_speeds_clear()
table = self.table_list_critical_speeds
critical_speeds = self.critital_speeds_calculate(modes_dict,ratio)
length = len(critical_speeds)
table.setRowCount(length)
for index in range(length):
table.setItem(index,0,QTableWidgetItem(str(critical_speeds[index])))
table.item(index,0).setFlags(QtCore.Qt.ItemIsEditable)
table.item(index,0).setTextAlignment(QtCore.Qt.AlignHCenter)
# BLOCK OF CREATING HANDLERS BUTTONS BEGIN!!
# create button calculate handler
def btn_calculate_handler(self):
self.terminal.clear()
os.system('cls')
self.btn_calculate.setEnabled(False)
self.toggle_visibility_widgets(False)
self.text_box_ratio_edit.setText('1')
self.plot_graph.clear()
QtTest.QTest.qWait(2000)
self.terminal.insertPlainText('Calculation start\n')
QtTest.QTest.qWait(2000)
self.terminal.insertPlainText('Clear Working Folder\n')
clear_working_folder(BASE_DIR,RESULT_FILENAME)
QtTest.QTest.qWait(3000)
self.terminal.insertPlainText('Clear Working Folder finished\n')
modes_dict = self.get_frequencies_from_vtu(self.calculate_damping_frequencies)
self.plot_campell_diagramm(modes_dict)
self.insert_critical_speeds_into_table(modes_dict)
self.btn_calculate.setEnabled(True)
self.terminal.insertPlainText('Calculation finished\n')
# create button replot_graph_and_recalculate_critical_speeds handler
def button_replot_graph_and_recalculate_critical_speeds_handler(self):
try:
self.terminal.insertPlainText(f'plot and calculate critical speeds with ratio {self.text_box_ratio_edit.text()} of critical line \n')
self.plot_graph.clear()
modes_dict = self.get_frequencies_from_vtu(self.empty_function)
self.plot_campell_diagramm(modes_dict,float(self.text_box_ratio_edit.text()))
self.insert_critical_speeds_into_table(modes_dict,float(self.text_box_ratio_edit.text()))
except ValueError:
self.message_box = QMessageBox(self)
self.message_box.setText('Invalid ratio')
self.message_box.exec_()
# BLOCK OF CREATING HANDLERS BUTTONS END!!
def application():
app = QApplication(sys.argv)
window = Window()
window.show()
sys.exit(app.exec_())
if __name__ == '__main__':
application()
Сохраните его с названием Window.py
Запустите на расчет файл Window.py. Выберите Run - Run Module.
После того как расчет прошел появится диалоговое окно для расчета:
Нажмите на Calculate fidesys script. В результате построится диаграмма Кемпбелла и посчитаются критические скорости для коэффициента критической линии 1:
В поле RATIO можно поменять коэффициент критической линии и диаграмма Кэмпбелла с табличкой расчета критических скоростей перерасчитается заново (с помощью кнопки Replot).
