Оптимизационная задача с использованием Fidesys Python API

Для проведения оптимизационного расчета необходимо, чтобы были выполнены следующие условия:

Установлен Python v.3.8 и выше;
Установлена библиотека vtk для Python;
Установлена библиотека numpy для Python.

Чтобы выполнить эти условия, необходимо сделать следующее:

Скачать Python 3.8 и выше с сайта python.org и установить.
Открыть командную строку Windows (cmd.exe) и прописать в ней:
pip3 install numpy (затем нажать Enter и дать установиться);
pip3 install vtk (затем нажать Enter и дать установиться).

После того как все подготовительные действия выполнены, можно приступать к решению задачи.

Рассматривается задача оптимизации диаметра основания столба рекламного щита, нагруженного ветровой нагрузкой.

Построение модели

1. Создайте параллелепипед.

На панели команд выберите модуль построения объёмной геометрии (Режим — Геометрия, Объект — Объём, Действие — Создать).

Из списка геометрических примитивов выберите Параллелепипед.

Задайте следующие параметры :

X (ширина): 20;
Y (высота): 0.5;
Z (глубина): 10.

Нажмите Применить.

2. Далее получившийся объем необходимо переместить.

На панели команд выберите модуль построения объёмной геометрии (Режим — Геометрия, Объект — Объём, Действие — Преобразовать).

Из списка операций выберите Переместить.

Задайте следующие параметры:

ID объема(ов): 1;
Выбрать метод: Расстояние;
Расстояние по Х: 0;
Расстояние по Y: 0;
Расстояние по Z: 30.

Нажмите Применить.

3. Создайте столб в виде усеченного конуса.

Из списка геометрических примитивов выберите Конус.

Задайте следующие параметры:

Высота: 25;
Верхний радиус: 0.25;
Круговой;
Радиус: 0.3.

Нажмите Применить.

4. Далее получившийся усеченный конус необходимо переместить.

Из списка операций выберите Переместить.

Задайте следующие параметры:

ID объема(ов): 2;
Выбрать метод: Расстояние;
Расстояние по Х: 0;
Расстояние по Y: 0;
Расстояние по Z: 12.5.

Нажмите Применить.

5. Создайте общие поверхности для построения совместной сетки.

На панели команд выберите модуль построения объёмной геометрии (Режим — Геометрия, Объект — Объём, Действие — Отпечатать/срастить).

Из списка операций выберите Отпечать и срастить.

Задайте следующие параметры:

ID объема(ов): all.

Нажмите Применить.

Построение сетки

1. Создайте тетраэдральную сетку.

На панели команд выберите модуль построения объёмной сетки (Режим — Сетка, Объект — Объемная, Действие — Построение сетки).

Из списка операций выберите Тетраэдральная.

Задайте следующие параметры:

Выбор объёмов: all.

Нажмите Построить сетку.

Если всё было сделано правильно, то вы должны увидеть такую модель:

Задание материала и свойств блока

1. Создайте материал.

На панели команд выберите модуль задания свойств материала (Режим — Материал, Объект — Управление материалами).

В открывшемся окне Управление материалами из третьей колонки во вторую перенесите материал «Углеродистая сталь».

Нажмите Применить.

Закройте окно Управление материалами.

2. Создайте блок.

На панели команд выберите модуль управления блоками (Режим — Блоки, Объект — Блок, Действие — Добавить сущность в блок).

Задайте следующие параметры:

Список сущностей: Объем;
ID объекта(ов): all.

Нажмите Применить.

3. Задайте параметры блока.

На панели команд выберите модуль управления блоками (Режим — Блоки, Объект — Блок, Действие — Свойства/параметры блока).

Задайте следующие параметры:

ID блока(ов): 1;
Категория: Объемное тело;
Материал: Углеродистая сталь;
Система координат: Глобальная декартова;
Порядок: 1.

Нажмите Применить.

Задание граничных условий

1. Закрепите по всем перемещениям основание столба.

На панели команд выберите модуль задания граничных условий (Режим — Граничные условия, Объект — Перемещение, Действие — Создать).

Задайте следующие параметры:

Список сущностей: Поверхность;
ID объекта(ов): 8;
Степени свободы: Закрепить все.

Нажмите Применить.

2. Задайте распределенную ветровую нагрузку на поверхности щита p = 230 Н/м².

На панели команд выберите модуль задания граничных условий (Режим — Граничные условия, Объект — Распределенная сила, Действие — Создать).

Задайте следующие параметры:

Список сущностей: Поверхность;
ID объекта(ов): 3;
Тип силы: Распределенная;
Сила: 230;
Вектор направления X: 0;
Вектор направления Y: 1;
Вектор направления Z: 0.

Нажмите Применить.

3. Добавьте действие гравитации.

На панели команд выберите модуль задания граничных условий (Режим — Граничные условия, Объект — Гравитация, Действие — Создать).

Задайте следующие параметры:

Список сущностей: Глобально;
Направления: Z;
Значение: -9.81.

Нажмите Применить.

Задание настроек решателя

1. Задайте настройки для статического решения задачи в 3D постановке.

На панели команд выберите модуль настройки расчёта (Режим — Настройки расчёта, Настройки расчета — Статический, Статический — Общие).

Задайте следующие параметры:

Размерность: 3D;
Модель: Упругость.

Нажмите Применить.

Нажмите Начать расчет.

2. В появившемся окне выберите директорию, в которой будет сохранён результат

3. В случае успешно проведённого расчёта в консоли отобразится сообщение: Calculation finished successfully at "date" "time".

Получение и преобразование скрипта модели

1. Извлеките скрипт модели из истории.

Перейдите в командную строку и переключите вкладку на "История", где вы увидите скрипт сформированной вами модели:

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
calculation start path 'D:/Fidesys/example_40.pvd'

Правой кнопкой мыши кликните в любое место в командной строке и выберите "Выделить всё", затем повторно кликните правой кнопкой по выделенному скрипту и выберите "Копировать".

Таким образом, вы скопировали скрипт в буфер обмена.

2. Преобразуйте скрипт в синтаксис Python.

Откройте редактор журнала и вставьте скопированный ранее скрипт в его окно.

Преобразуйте скрипт в синтаксис Python через меню Инструменты - Интерпретатор - Python.

Если всё сделано правильно, то вы получите следующий скрипт в окне:

Скопируйте из редактора журнала получившийся скрипт.

Создание и запуск скрипта Python

1. Создайте файл скрипта Python.

Запустите Python IDLE, в меню выберите File - New File, после чего откроется окно для редактирования скрипта.

2. Скопируйте приведенный ниже скрипт в открывшееся пустое окно.

Данный скрипт Python уже содержит ту часть скрипта модели Fidesys, которую мы получили ранее. Место вставки скрипта модели Fidesys помечено соответствующими комментариями.

Обратите внимание, что варьирование диаметра основания столба происходит путем внесений изменений в команду создания конуса:

- начальный вид команды: fidesys.cmd("create frustum height 25 radius 0.25 top 0.25");

- вид изменяемой команды: fidesys.cmd("create frustum height 25 radius "+str(r)+"top 0.25").

Вставка "+str(r)+" добавляет в разрыв текстовой команды значение радиуса.

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-9.0'       # Расположение Фидесиса
base_dir = os.path.dirname(os.path.abspath(__file__))            # Директория где лежит скрипт    
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
        
r = 0.25             # Начальный радиус основания
isOptimized = False  # Изначально False - начальная конструкция не оптимизирована
iteration = 1        # Начальное значение счетчика итераций
limit = 100          # Ограничение на число итераций
        
while isOptimized == False and iteration <= limit:
    print("Итерация № ",iteration) # Пишем в консоль какой проход
    print("Диаметр ", 2*r)
    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("Начинаем расчет в файл " + output_pvd_path)               # Выводим в консоль директорию и файл сохранения
    fidesys.cmd("calculation start path '" + output_pvd_path + "'")  # Просим Фидесис начать расчет в указанную директорию
        	
    print("                  ")
    print("Расчет завершен!\n")
    print("Начинается чтение результатов.")
        
    reader = vtk.vtkXMLUnstructuredGridReader()                             #Подключаем читалку
    filename = os.path.join(str(base_dir)+r"\1\case1_step0001_substep0001.vtu") # Указываем путь к файлу
    print("Читаем результаты из ",filename)                                 # Пишет откуда берем результаты
    reader.SetFileName(filename)                                            # Подключаем путь к читалке и читаем
    reader.Update()                                                         # Needed because of GetScalarRange
    grid = reader.GetOutput()                                               # Забираем выходные данные
    point_data = grid.GetPointData()                                        # Забираем данные для точек
    print("Чтение результатов произведено.\n")
        
    arrayOfStress = vtk_to_numpy(point_data.GetArray("Stress")) # Считываем напряжения из массива результатов
    node_id = vtk_to_numpy(point_data.GetArray("Node ID"))      # Считываем номера узлов из массива результатов
        
    print("Начинаем поиск перенапряженных узлов")
    for i in range(len(arrayOfStress)):         # range(len(array...)) создает массив чисел для счетчика цикла for
            if arrayOfStress[i][6] > 106e6:     # Проверяем напряжения по Мизесу в узлах
                overstressed.append(node_id[i]) # Заполняем массив номерами перенапряженных узлов
        
    if len(overstressed) == 0: # Проверяется размер массива перенапряженных узлов, если он 0 то
        isOptimized = True     # приравниваем переменную isOptimized=True, чтобы выйти из цикла
        print("Конструкция оптимизирована!")
    else:
        print("Перенапряженных узлов: ",len(overstressed)) # Выводим информацию о количестве перенапряженных узлов
        print("                  ")
        r = r + 0.025    # Увеличиваем радиус на 0.05
        iteration = iteration + 1 # Увеличиваем значение счетчика проходов            		
fc.delete_application() # Команда удаления задачи из памяти для версии 5.2 (для 5.0 и ниже заменить на fc.deleteApplication())
print("                  ")
if isOptimized == True:
    print("Готово! Оптимальный диаметр не менее ", 2*r)
else:
    print("Расчет остановлен, т.к. достигнуто предельное число итераций ", limit)

3. Запустите скрипт.

В меню выберите Run - Run Module. Когда система попросит сохранить этот файл, сохраните его в папку Example, созданную в директории, в пути которой нет кириллических символов, в избежание ошибок, и назовите его optimization.

В процессе и по завершению расчета вы будете наблюдать следующие результаты в консоли.

Сохранение результатов произойдет в ту папку, в которой лежит файл py-скрипта. По завершению расчета вы сможете открыть и посмотреть файл результатов 1.pvd.