Моделирование роста трещины при хрупком разрушении в 2D постановке в Fidesys 5.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-5.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_paths = [os.path.join(fidesys_path, 'preprocessor', 'bin'),
os.path.join(fidesys_path, 'preprocessor', 'bin', 'plugins'),
os.path.join(fidesys_path, 'preprocessor', 'bin', 'acis', 'code', 'bin'),
os.path.join(fidesys_path, 'preprocessor', 'structure')] # Директории, где препроцессор
calc_path = os.path.join(fidesys_path,'bin', 'FidesysCalc.exe') # Директория, где решатель
for path in prep_paths:
os.environ['PATH'] += os.pathsep + path # Добавление пути к препроцессору в PATH
os.add_dll_directory(path)
sys.path.append(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_paths[0]) # !!!!!Инициализация для версий 5.1 (для 5.0 и ниже заменить на fc.initApplication(prep_paths[0]) )
fc.start_up_no_args() # Запуск обязательного компонента Фидесис fc
deleteNd=[] #Пустой массив узлов для удаления
deleteEl=[] #Пустой массив элементов для удаления
for epoch in range(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('create surface rectangle width 1 height 1 zplane ')
fidesys.cmd('create surface circle radius 0.05 zplane ')
fidesys.cmd('move Surface 2 y -0.3 include_merged preview ')
fidesys.cmd('move Surface 2 y -0.3 include_merged ')
fidesys.cmd('subtract surface 2 from surface 1 ')
fidesys.cmd('create vertex -0.02312 0.5 0 on curve 1 ')
fidesys.cmd('create vertex 0.04413 0.5 0 on curve 1 ')
fidesys.cmd('create vertex -0.002101 -0.2 0 on surface 3 ')
fidesys.cmd('create surface vertex 7 9 8 ')
fidesys.cmd('subtract surface 4 from surface 3')
fidesys.cmd('webcut body 1 with plane xplane offset 0.2 imprint merge ')
fidesys.cmd('webcut body 4 with plane xplane offset -0.2 imprint merge ')
fidesys.cmd('surface 9 6 size auto factor 1')
fidesys.cmd('mesh surface 9 6 ')
fidesys.cmd('curve 11 interval 60')
fidesys.cmd('curve 11 scheme bias factor 1.05')
fidesys.cmd('curve 10 interval 60')
fidesys.cmd('curve 10 scheme bias factor 0.95')
fidesys.cmd('curve 6 interval 60')
fidesys.cmd('curve 6 scheme equal')
fidesys.cmd('surface 8 scheme trimesh geometry approximation angle 15 ')
fidesys.cmd('Trimesher surface gradation 1.3')
fidesys.cmd('Trimesher geometry sizing on')
fidesys.cmd('surface 8 size 0.005')
fidesys.cmd('mesh surface 8')
fidesys.cmd('renumber Node all start_id 1 uniqueids')
fidesys.cmd('renumber element all start_id 1 uniqueids ')
fidesys.cmd('create material 1 from \'Углеродистая сталь\'')
fidesys.cmd('set duplicate block elements off')
fidesys.cmd('block 1 add surface all')
fidesys.cmd('block \'Block 1\' material 1 cs 1 element plane order 1')
fidesys.cmd('create displacement on curve 21 dof all fix ')
fidesys.cmd('create pressure on curve 2 4 magnitude -1000000')
fidesys.cmd('analysis type static elasticity dim2 planestress')
# ---------Конец вставляемого скрипта из Фидесис-------------
print("Модель создана")
if epoch != 0:
# здесь мы удаляем элементы, которые выбрали в предыдущем шаге расчета
print("Начинаем удаление элементов")
fidesys.cmd('set dev on')
for element in deleteEl:
#print("delete Tet " + str(tetra))
fidesys.cmd("delete Element " + str(element))
# сообщаем в кубит-файл что здесь мы пытаемся начать расчет с данными предыдущего шага
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()
element_q = cubit.get_element_count()
print("it includes","nodes={},elements={}".format(nodes_q,element_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_step0001_substep0001.vtu'.format(epoch)) # Пишет откуда берем результаты
filename = os.path.join(base_dir, str(epoch+1), 'case1_step0001_substep0001.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] > 100e6: # Проверяем 1 главные напряжения в узлах (если хотите по Мизесу, то arrayOfStress[point][6])
deleteNd.append(node_id[point]) # Заполняем массив номерами перенапряженных узлов
print("Узлов к удалению: ",len(deleteNd))
watcher = -1 #следилка за номером элемента, чтобы не дублировать запись об удалении
for element in range(1,element_q+1):
nodes = cubit.get_expanded_connectivity("Element", element)
for node in nodes:
for ndToDelete in deleteNd:
if (node == ndToDelete):
if (watcher != element):
#print("Element ", element, "will be deleted")
watcher = element
deleteEl.append(element)
print("Элементов к удалению ",len(deleteEl))
print(" ")
print("--------------------------------------")
fc.delete_application() #Команда удаления задачи из памяти для версии 5.1 (для 5.0 и ниже заменить на 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_step0001_substep0001.vtu" name="Name: Part: Material:ELASTIC"/>\n')
f.write(' </Collection>\n')
f.write('</VTKFile>\n')
print("done!")
fidesys