Расчет разрушения контейнера НЗК при ударе об землю в Fidesys 4.1 с использованием Python API

Рассматривается задача падения контейнера для хранения РАО на землю или иное жесткое основание.

Что минимально необходимо для использования:

  • 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 и дать установиться)

Далее необходимо создать модель и подготовить ее для расчета.

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

Отсеките один из углов.

Отсеките оставшиеся три угла.

Удалите отсеченный материал.

Создайте уменьшенную копию внутри.

Переместите внутренний объем.

Удалите «перекрытие».

Создайте объем для вычитания.

Определите номера поверхностей.

Переместите объем через перемещение поверхностей.

Вычтите нижний объем.

Отсеките края.

Отсеките углы.

Отсеките часть объема 16.

Отсеките части остальных объемов.

Удалите внутренние объемы.

Объедините все части контейнера.

Поверните контейнер на 45 градусов вокруг его оси.

Наклоните контейнер.

Создайте объем для контакта.

Переместите контактный объем.

Уменьшите размер задачи вдвое через отсечение половины модели.

Удалите лишние объемы.

Создайте общие плоскости РАО и НЗК для совместности сетки.

Для удобства «Подрежем ногу».

Удалите ненужную часть.

Задайте объемам метод сетки.

Постройте сетку сначала на РАО.

Затем постройте сетку на контейнере и поверхности.

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

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

Создайте 1й блок - с контейнером и РАО.

Создайте 2й блок - с объемом, который моделирует основание.

Присвойте свойства блокам.

Добавьте закрепление нижней плоскости основания.

Добавьте ГУ симметрии.

Создайте контактную пару.

Задайте нагрузку в виде гравитации.

Задайте начальную скорость падения.

Сблизьте контактные поверхности.

Задайте настройки расчета и запустите на расчет.

Выберите место сохранения и дождитесь окончания расчета, чтобы понять, что модель работает и контакт происходит.

Посмотрите результаты в постпроцессоре на предмет корректности происходящего.

Если все происходит правильно, то вернитесь в препроцессор и откройте снизу вкладку "История".

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

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

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

Должен получиться такой набор команд.

Далее необходимо создать новый файл Python.

И скопировать в открывшееся пустое окно следующий скрипт для версий 4.1

# Импортируем нужные библиотеки
import vtk # Библиотека работы с выходными данными
from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk # Модуль для преобразования результатов
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') # Файл ссылок на результаты
result_dir = os.path.join(base_dir, 'results') # Директория с результатами
prep_path = os.path.join(fidesys_path, 'preprocessor', 'bin') # Директория, где препроцессор
calc_path = os.path.join(fidesys_path,'bin', 'FidesysCalc.exe') # Директория, где решатель
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

deleteNd=[] #Пустой массив узлов для удаления
deleteEl=[] #Пустой массив элементов для удаления

#Функция, определяющая сколько добавить нулей к имени выходного файла т.к. их число зависит от количества цифр в номере шага
def zerosNumber(epoch):
    n = len(str(epoch))
    st = ""
    for i in range(n):
        st=st+str(0)
    return st

for epoch in range(0, 12): # Задаем количество шагов
    
    print("starting epoch "+str(epoch))
    fidesys.cmd('#--------------------------------------------------------')
    fidesys.cmd('#note for cubit.jou: prepairing for epoch '+str(epoch)+'(step '+str(epoch+1)+')')
######################## Начиная отсюда вставляем код модели ########################
    fidesys.cmd('reset')
    fidesys.cmd('brick x 1.640 y 1.640 z 1.355')
    fidesys.cmd('webcut volume all with plane xplane offset 1 rotate 45 about z center 0 0 0  ')
    fidesys.cmd('webcut volume all with plane xplane offset -1 rotate 45 about z center 0 0 0  ')
    fidesys.cmd('webcut volume all with plane xplane offset -1 rotate -45 about z center 0 0 0  ')
    fidesys.cmd('webcut volume all with plane xplane offset 1 rotate -45 about z center 0 0 0  ')
    fidesys.cmd('delete volume 1 2 4 3  ')
    fidesys.cmd('Volume 5  copy scale x 0.8 y 0.8 z 0.6 ')
    fidesys.cmd('move Volume 6 z 0.05 include_merged ')
    fidesys.cmd('remove overlap volume 6 5 modify larger')
    fidesys.cmd('Volume 6 copy scale x 1 y 1 z 0.2 ')
    fidesys.cmd('move Surface 79 location surface 44 except x y include_merged ')
    fidesys.cmd('subtract volume 9 from volume 8  ')
    fidesys.cmd('webcut volume 10  with plane from surface 92  ')
    fidesys.cmd('webcut volume 10  with plane from surface 87  ')
    fidesys.cmd('webcut volume 11  with plane from surface 87  ')
    fidesys.cmd('webcut volume 11  with plane from surface 93  ')
    fidesys.cmd('webcut volume 11  with plane from surface 89  ')
    fidesys.cmd('webcut volume 10 12 13 11  with plane from surface 90  ')
    fidesys.cmd('webcut volume 18  with plane xplane offset 1,25 ')
    fidesys.cmd('webcut volume 16  with plane xplane offset 1.25 rotate -45 about z ')
    fidesys.cmd('webcut volume 18  with plane xplane offset -1.25 rotate -45 about z ')
    fidesys.cmd('webcut volume 17  with plane xplane offset 1.25 rotate 45 about z ')
    fidesys.cmd('webcut volume 15  with plane xplane offset -1.25 rotate 45 about z ')
    fidesys.cmd('delete volume 15 19 21 18  ')
    fidesys.cmd('unite volume 16 12 14 10 22 11 20 13 17  ')
    fidesys.cmd('imprint volume all ')
    fidesys.cmd('merge volume all ')
    fidesys.cmd('rotate Volume all angle 45  about Z include_merged ')
    fidesys.cmd('rotate Volume all angle 55  about Y include_merged ')
    fidesys.cmd('brick x 0.5 y 0.5 z 0.1')
    fidesys.cmd('move Surface 229 midpoint location curve 69  except x y include_merged ')
    fidesys.cmd('webcut volume all with plane yplane ')
    fidesys.cmd('delete volume 25 26 27 ')
    fidesys.cmd('imprint volume 6 23  ')
    fidesys.cmd('merge volume 6 23 ')
    fidesys.cmd('webcut volume 23  with plane zplane offset -1.2 ')
    fidesys.cmd('delete volume 28  ')
    fidesys.cmd('volume 6 23 24 scheme tetmesh')
    fidesys.cmd('volume 6  size auto factor 5')
    fidesys.cmd('mesh volume 6 ')
    fidesys.cmd('volume 23  size auto factor 4')
    fidesys.cmd('mesh volume 23 ')
    fidesys.cmd('volume 24  size auto factor 8')
    fidesys.cmd('mesh volume 24 ')
    fidesys.cmd('create material 1')
    fidesys.cmd('modify material 1 name \'КонтейнерРАО\'')
    fidesys.cmd('modify material 1 set property \'DENSITY\' value 2400')
    fidesys.cmd('modify material 1 set property \'POISSON\' value 0.2')
    fidesys.cmd('modify material 1 set property \'MODULUS\' value 3e+10')
    fidesys.cmd('create material 2')
    fidesys.cmd('modify material 2 name \'Поверхность\'')
    fidesys.cmd('modify material 2 set property \'POISSON\' value 0')
    fidesys.cmd('modify material 2 set property \'DENSITY\' value 1e-10')
    fidesys.cmd('modify material 2 set property \'MODULUS\' value 1e+30')
    fidesys.cmd('set duplicate block elements off')
    fidesys.cmd('block 1 add volume 6 23 ')
    fidesys.cmd('set duplicate block elements off')
    fidesys.cmd('block 2 add volume 24 ')
    fidesys.cmd('block \'Block 1\' material 1 cs 1 element solid order 1')
    fidesys.cmd('block \'Block 2\' material 2 cs 1 element solid order 1')
    fidesys.cmd('create displacement  on surface 274  dof all fix  ')
    fidesys.cmd('create displacement  on surface 271 284 235  dof 2 fix  ')
    fidesys.cmd('create contact master surface 272 slave surface 283 285 281 tolerance 0.0005 type general friction 0.0 offset 0.0 ignore_overlap off method auto')
    fidesys.cmd('create gravity global')
    fidesys.cmd('modify gravity 1 dof 3 value -9.81')
    fidesys.cmd('create initial velocity on volume 6 23 ')
    fidesys.cmd('modify initial velocity 1 dof 3 value -1')
    fidesys.cmd('move Surface 281  location surface 272  except x y include_merged ')
    fidesys.cmd('analysis type dynamic elasticity dim3 preload off')
################### До сюда вставляется код модели, а на след.строке настройки одного шага ######################
    fidesys.cmd('dynamic method full_solution scheme implicit steps 1 newmark_gamma 0.005 maxtime 0.001')
    print("Модель создана")
###################
    if epoch != 0:
        # Считывание файла результатов предыдущего шага
        input_vtu_path = os.path.join(base_dir, str(epoch), 'case1_step'+zerosNumber(epoch)+'001.vtu'.format(epoch))
        if not os.path.exists(input_vtu_path): # Проверка наличия директории
            print('Not vtu {}'.format(input_vtu_path))
            sys.exit(1)
            
        # Направляем команду импорта предыдущего шага
        import_vtu = 'import vtu "{}" mesh_geometry linear result_as_initial keep_data skip_mesh'.format(input_vtu_path)
        print(import_vtu)
        fidesys.cmd(import_vtu)
        
        # Направляем команду для следующего шага в динамическом расчете
        settings = "dynamic method full_solution scheme implicit maxtime " + \
              str(0.001+epoch*0.001) + " steps "+str(epoch+1) + \
              " newmark_gamma 0.005 inittime "+str(epoch*0.001)+" initstep "+ str(epoch)
        print(settings)
        fidesys.cmd(settings)

        # здесь мы удаляем элементы, которые выбрали в предыдущем шаге расчета
        print("Начинаем удаление элементов")
        for tetra in deleteEl:
            fidesys.cmd('set dev on')
            #print("delete Tet " + str(tetra))
            fidesys.cmd("delete Tet " + str(tetra))

        # сообщаем в кубит-файл что здесь мы пытаемся начать расчет с данными предыдущего шага
        fidesys.cmd('#note for cubit.jou: here we trying to start solution using data of step '+str(epoch)+'(epoch '+str(epoch-1)+')')
    else:
        fidesys.cmd('#note for cubit.jou: here we trying to start solution of 0 epoch (1-st step to 1.pvd)')

    # Смотрим сколько узлов и элементов в модели
    nodes_q = cubit.get_node_count()
    tets_q = cubit.get_tet_count()
    print("it includes","nodes={},tets={}".format(nodes_q,tets_q))
    print("                  ")
    
    # Задаем путь - куда совершать расчет
    output_pvd_path = os.path.join(base_dir + "\\" + str(epoch+1)+".pvd")
    print("strarting calculation to " + output_pvd_path)

    #Направляем команду для совершения расчета в указанный путь
    fidesys.cmd("calculation start path '" + output_pvd_path + "'")
    print("                  ")

    print("Расчет удачно завершен!")
    print("                  ")

    reader = vtk.vtkXMLUnstructuredGridReader() #Подключаем читалку
    print("Читаем результаты из "+base_dir+"\\"+str(epoch+1)+"\\"+'case1_step'+zerosNumber(epoch+1)+'001.vtu'.format(epoch)) # Пишет откуда берем результаты
    filename = os.path.join(base_dir, str(epoch+1), 'case1_step'+zerosNumber(epoch+1)+'001.vtu'.format(epoch)) # Указываем путь к файлу
    reader.SetFileName(filename) # Подключаем путь к читалке и читаем
    print("Результаты прочитаны!")
    print("                  ")
    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][7] > 30e6: # Проверяем 1 главные напряжения в узлах (если хотите по Мизесу, то arrayOfStress[point][6])
                        deleteNd.append(node_id[point]) # Заполняем массив номерами перенапряженных узлов
    print("Узлов к удалению: ",len(deleteNd))

    
    watcher = -1 #следилка за номером элемента, чтобы не дублировать запись об удалении
    for tetra in range(1,tets_q+1):
                nodes = cubit.get_expanded_connectivity("Tet", tetra)
                for node in nodes: 
                    for ndToDelete in deleteNd: 
                        if (node == ndToDelete):
                            if (watcher != tetra):
                                #print("Tet ", tetra, "will be deleted")
                                watcher = tetra
                                deleteEl.append(tetra)
    print("elements to be deleted ",len(deleteEl))
    print("                  ")
    print("--------------------------------------")

fc.deleteApplication() #Удаляем из памяти выполненную задачу
print("Готово!")

print("--------------------------------------")
print("Делаем файл ссылок на все результаты - all-results.pvd")
with open(str(base_dir)+ "\\all-results.pvd", "w") as f:
   
# Пишем в файл ссылки на результаты
    f.write('<?xml version="1.0"?>\n')
    f.write('<VTKFile type="Collection" version="0.1" byte_order="LittleEndian" compressor="vtkZLibDataCompressor">\n')
    f.write('  <Collection>\n')
    f.write('    <DataSet timestep="0.0" part="0" file="1/case1_step0000.vtu" name="Name: Part: Material:ELASTIC"/>\n')
    for i in range(1, epoch+2):
        f.write('    <DataSet timestep="'+str(i*0.001)+'" part="0" file="'+str(i)+'/case1_step'+zerosNumber(i)+'001.vtu" name="Name: Part: Material:ELASTIC"/>\n')
    f.write('  </Collection>\n')
    f.write('</VTKFile>\n')
print("done!")

Запускаем на расчет.

Если все сделано правильно, то вы увидите такой процесс на экране.

По завершении расчета вы увидите следующее.

После чего вы сможете посмотреть результаты в all-results.pvd.