Рассматривается задача оптимизации диаметра основания столба рекламного щита, нагруженного ветровой нагрузкой.
Свойства материала (для простоты) берутся из библиотеки Фидесис следующими:
Что минимально необходимо для использования:
- Python v.3.8
- Библиотека vtk – для чтения результатов
- Библиотека numpy – для преобразования результатов
- Библиотеки sys и os – для работы с путями к файлам
- Библиотеки cubit и fidesys – для препроцессинга данных
Подготовка:
- Скачать Python 3.8 с сайта python.org и установить.
- Скачать Fidesys 4.1 с сайта cae-fidesys.com и установить.
- Открыть командную строку Windows (cmd.exe) и прописать в ней:
- pip3 install vtk (затем нажать Enter и дать установиться)
- pip3 install numpy (затем нажать Enter и дать установиться)
Далее необходимо создать модель и подготовить ее для расчета.
Создаем сам щит.
Перемещаем его.
Создаем столб в виде усеченного конуса.
Перемещаем его.
Создаем общие поверхности для того, чтобы сетка совместилась.
Строим простую тетраэдральную сетку.
Если всё было сделано правильно, до вы должны увидеть такую модель:
Создаем материал путем добавления материала из базы данных.
Создаем блок из всех объемов, чтобы передать им свойства модели.
Присваиваем свойства блоку.
Создаем граничные условия основанию столба.
Задаем ветровую нагрузку
Добавляем действие гравитации.
Задаем настройки расчета.
В итоге получится следующий скрипт, который можно будет увидеть во вкладке "История" командной строки:
reset
brick x 20 y 0.5 z 10
move Volume 1 x 0 y 0 z 30 include_merged
create frustum height 25 radius 0.25 top 0.25
move Volume 2 x 0 y 0 z 12.5 include_merged
undo group begin
imprint volume all
merge volume all
undo group end
volume all scheme tetmesh
mesh volume all
create material 1 from 'Углеродистая сталь'
set duplicate block elements off
block 1 add volume all
block 'Block 1' material 1 cs 1 element solid order 1
create displacement on surface 8 dof all fix
create distributed force on surface 3 force value 230 moment value 0 direction 0 1 0 specific
create gravity global
modify gravity 1 dof 3 value -9.81
analysis type static elasticity dim3
Затем его необходимо скопировать для преобразования в синтаксис Питона.
Далее открываем редактор журнала и вставляем в пустое окно текст скрипта.
Затем преобразуем в синтаксис Питона.
Получим следующий скрипт:
Далее необходимо создать новый файл Python.
И скопировать в открывшееся пустое окно следующий скрипт:
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-4.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
import cubit # Библиотека препроцессинга
import fidesys # Библиотека Фидесиса
cubit.init([""]) # Инициализация препроцессора
fc = fidesys.FidesysComponent() # Создание обязательного компонента Фидесис fc
fc.initApplication(prep_path)
fc.start_up_no_args() # Запуск обязательного компонента Фидесис fc
r = 0.25 # Начальный радиус основания
print("Начальный диаметра ", 2*r)
isOptimized = False # Изначально False - начальная конструкция не оптимизирована
prohod = 1 # Начальное значение счетчика проходов(итераций)
while isOptimized == False:
print("Проход № ",prohod) # Пишем в консоль какой проход
overstressed = [] # Создаем пустой массив для заполнения перенапряженными узлами
#--------- Сюда вставляется геометрический скрипт с настройками расчета ----------
output_pvd_path = os.path.join(base_dir + "\\" + "1.pvd")
print("strarting calculation to " + output_pvd_path)
fidesys.cmd("calculation start path '" + output_pvd_path + "'")
print(" ")
print("Расчет успешно завершен!")
print(" ")
reader = vtk.vtkXMLUnstructuredGridReader() #Подключаем читалку
print("Читаем результаты из ",str(base_dir)+r"\1\case1_step01_substep0001.vtu") # Пишет откуда берем результаты
filename = os.path.join(str(base_dir)+r"\1\case1_step01_substep0001.vtu") # Указываем путь к файлу
reader.SetFileName(filename) # Подключаем путь к читалке и читаем
reader.Update() # Needed because of GetScalarRange
grid = reader.GetOutput() # Забираем выходные данные
point_data = grid.GetPointData() # Забираем данные для точек
arrayOfStress = vtk_to_numpy(point_data.GetArray("Stress")) # Считываем напряжения из массива результатов
node_id = vtk_to_numpy(point_data.GetArray("Node ID")) # Считываем номера узлов из массива результатов
print("Начинаем поиск перенапряженных узлов")
print(" ")
for point in range(len(arrayOfStress)):
if arrayOfStress[point][6] > 106e6: # Проверяем напряжения по Мизесу в узлах
overstressed.append(node_id[point]) # Заполняем массив номерами перенапряженных узлов
if len(overstressed) == 0: # Проверяется размер массива перенапряженных узлов, если он 0 то
isOptimized = True # приравниваем переменную isOptimized=True, чтобы выйти из цикла
print("Конструкция оптимизирована!")
else:
print("Перенапряженных узлов: ",len(overstressed)) # Выводим информацию о количестве перенапряженных узлов
print(" ")
r = r + 0.05 # Увеличиваем радиус на 0.05
prohod = prohod + 1 # Увеличиваем значение счетчика проходов
fc.deleteApplication() #Удаляем из памяти выполненную задачу
print(" ")
print("Готово! Оптимальный диаметр не менее ", 2*r)
Затем необходимо вставить вместо текста#--------- Сюда вставляется геометрический скрипт с настройками расчета ----------непосредственно тот преобразованный текст модели, который мы получили ранее, затем делаем вставку "+str(r)+" вместо первого 0.25 в строчку fidesys.cmd("create frustum height 25 radius "+str(r)+"top 0.25"), чтобы у нас перебиралось значение нижнего радиуса конуса, после чего в итоге получится следующий скрипт:
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-4.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
import cubit # Библиотека препроцессинга
import fidesys # Библиотека Фидесиса
cubit.init([""]) # Инициализация препроцессора
fc = fidesys.FidesysComponent() # Создание обязательного компонента Фидесис fc
fc.initApplication(prep_path)
fc.start_up_no_args() # Запуск обязательного компонента Фидесис fc
r = 0.25 # Начальный радиус основания
print("Начальный диаметра ", 2*r)
isOptimized = False # Изначально False - начальная конструкция не оптимизирована
prohod = 1 # Начальное значение счетчика проходов(итераций)
while isOptimized == False:
print("Проход № ",prohod) # Пишем в консоль какой проход
overstressed = [] # Создаем пустой массив для заполнения перенапряженными узлами
# ---------Начало вставляемого скрипта из Фидесис-------------
fidesys.cmd("reset")
fidesys.cmd("brick x 20 y 0.5 z 10")
fidesys.cmd("move Volume 1 x 0 y 0 z 30 include_merged ")
fidesys.cmd("create frustum height 25 radius "+str(r)+"top 0.25")
fidesys.cmd("move Volume 2 x 0 y 0 z 12.5 include_merged ")
fidesys.cmd("undo group begin")
fidesys.cmd("imprint volume all ")
fidesys.cmd("merge volume all ")
fidesys.cmd("undo group end")
fidesys.cmd("volume all scheme tetmesh")
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("create displacement on surface 8 dof all fix ")
fidesys.cmd("create distributed force on surface 3 force value 230 moment value 0 direction 0 1 0 specific")
fidesys.cmd("create gravity global")
fidesys.cmd("modify gravity 1 dof 3 value -9.81")
fidesys.cmd("analysis type static elasticity dim3")
# ---------Конец вставляемого скрипта из Фидесис-------------
output_pvd_path = os.path.join(base_dir + "\\" + "1.pvd") # Объявляем директорию и файл сохранения
print("strarting calculation to " + output_pvd_path) # Выводим в консоль директорию и файл сохранения
fidesys.cmd("calculation start path '" + output_pvd_path + "'") # Просим Фидесис начать расчет в указанную директорию
print(" ")
print("Расчет успешно завершен!")
print(" ")
reader = vtk.vtkXMLUnstructuredGridReader() #Подключаем читалку
print("Читаем результаты из ",str(base_dir)+r"\1\case1_step0001_substep0001.vtu") # Пишет откуда берем результаты
filename = os.path.join(str(base_dir)+r"\1\case1_step0001_substep0001.vtu") # Указываем путь к файлу
reader.SetFileName(filename) # Подключаем путь к читалке и читаем
reader.Update() # Needed because of GetScalarRange
grid = reader.GetOutput() # Забираем выходные данные
point_data = grid.GetPointData() # Забираем данные для точек
arrayOfStress = vtk_to_numpy(point_data.GetArray("Stress")) # Считываем напряжения из массива результатов
node_id = vtk_to_numpy(point_data.GetArray("Node ID")) # Считываем номера узлов из массива результатов
print("Начинаем поиск перенапряженных узлов")
print(" ")
for point in range(len(arrayOfStress)):
if arrayOfStress[point][6] > 106e6: # Проверяем напряжения по Мизесу в узлах
overstressed.append(node_id[point]) # Заполняем массив номерами перенапряженных узлов
if len(overstressed) == 0: # Проверяется размер массива перенапряженных узлов, если он 0 то
isOptimized = True # приравниваем переменную isOptimized=True, чтобы выйти из цикла
print("Конструкция оптимизирована!")
else:
print("Перенапряженных узлов: ",len(overstressed)) # Выводим информацию о количестве перенапряженных узлов
print(" ")
r = r + 0.05 # Увеличиваем радиус на 0.05
prohod = prohod + 1 # Увеличиваем значение счетчика проходов
fc.deleteApplication() #Удаляем из памяти выполненную задачу
print(" ")
print("Готово! Оптимальный диаметр не менее ", 2*r)
Далее скрипт необходимо запустить и, когда система попросит, сохранить этот файл к примеру в папку Example, созданную на рабочем столе.
Если возникнут ошибки чтения/записи, то проверьте нет ли в тексте пути кириллических символов. Во избежание ошибок желательно, чтобы путь содержал только латинские буквы.
Запускаем.
После того как расчет прошел вы увидите:
И в папке с файлом скрипта будут лежать следующие файлы:
.svg)
.jpg)
.jpg)
