Определение долговечности алюминиевого образца при воздействии случайной широкополосной вибрации

Данная задача демонстрирует возможность описания усталостных свойств с помощью кривой зависимости напряжений от количества циклов до разрушения (кривая Велера, кривая S-N).

В качестве примера используется растянутый образец из сплава алюминия 6061-T6, подвергающийся воздействию широкополосной вибрации в течении 24 часов. Необходимо определить срок службы и коэффициент накопленной поврежденности.

Установите утилиту для создания виртуального окружения в python virtualenv:

pip install  virtualenv

Создайте произвольную папку (например, “fatique analysis”), создайте в ней два файла: calculate.py и fatique.py. Также там необходимо создать таблицу csv со спектром ускорения:

calculate.py:

import subprocess, os # импорт необходимых библиотек
FIDESYS_PATH = 'C:\Program Files\Fidesys\CAE-Fidesys-6.0' # путь к установке фидесиса
PVPYTHON_PATH  = os.path.join(FIDESYS_PATH,'postprocessor','pvpython.exe') # путь к pvpython
subprocess.run([PVPYTHON_PATH,'fatique.py'],check=True) # запуск файла fatique.py через консоль

fatique.py:

# блок для запуска виртуального окружения через скрипт питона
import os
BASE_DIR = os.getcwd() # Директория, где находится скрипт
activate_this = os.path.join(BASE_DIR,"venv","Scripts","activate_this.py")
with open(activate_this) as f:
    exec(f.read(), {'__file__': activate_this})
f.close()
# блок импорта необходимых библиотек
import meshio 
import sys 
import math
import numpy as np 
from paraview.simple import OpenDataFile, servermanager, RandomVibration, _select, ExtractSelection
from paraview.numpy_support import vtk_to_numpy 

# блок подключения библиотек предпроцессора фидесиса
FIDESYS_PREP_DIR = 'C:/Program Files/Fidesys/CAE-Fidesys-6.0/preprocessor/bin' 
prep_paths = [
FIDESYS_PREP_DIR,
os.path.join(FIDESYS_PREP_DIR, 'plugins'),
os.path.join(FIDESYS_PREP_DIR, 'acis', 'code', 'bin')
] # директории, где препроцессор
calc_path = os.path.join(FIDESYS_PREP_DIR, '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, fidesys
except ModuleNotFoundError:
    print(f'Модули fydesys, cubit не найдены в директории {FIDESYS_PREP_DIR}')  
    sys.exit(1) 

# блок дополнительных констант
DURATION_TIME = 86400 # Время длительности испытания на усталость (24 часа)
K = -0.102 # константа K для  log(sigma) = K*log(N) + B 
B = 8.75    # константа B для  log(sigma) = K*log(N) + B 
# Блок для работы с предпроцессором
cubit.init([""]) # инициализация функций препроцессора
fc = fidesys.FidesysComponent() # cоздание обязательного компонента Фидесис fc
fc.init_application(FIDESYS_PREP_DIR) # инициализация
fc.start_up_no_args() # запуск обязательного компонента Фидесис fc
fidesys.cmd('reset')
fidesys.cmd('create vertex 0 0 0')
fidesys.cmd('create vertex 0.03 0 0')
fidesys.cmd('create vertex 0.04 -0.0035 0')
fidesys.cmd('create vertex 0.12 -0.0035 0')
fidesys.cmd('create vertex 0.13 0 0')
fidesys.cmd('create vertex 0.16 0 0')
fidesys.cmd('create vertex 0 -0.02 0')
fidesys.cmd('create vertex 0.03 -0.02 0')
fidesys.cmd('create vertex 0.04 -0.0165 0')
fidesys.cmd('create vertex 0.12 -0.01655 0')
fidesys.cmd('create vertex 0.13 -0.02 0')
fidesys.cmd('create vertex 0.16 -0.02 0')
fidesys.cmd('create curve vertex 1 2')
fidesys.cmd('create curve vertex 2 3')
fidesys.cmd('create curve vertex 3 4')
fidesys.cmd('create curve vertex 4 5')
fidesys.cmd('create curve vertex 5 6')
fidesys.cmd('create curve vertex 6 12')
fidesys.cmd('create curve vertex 12 11')
fidesys.cmd('create curve vertex 11 10')
fidesys.cmd('create curve vertex 10 9')
fidesys.cmd('create curve vertex 9 8')
fidesys.cmd('create curve vertex 8 7')
fidesys.cmd('create curve vertex 7 1')
fidesys.cmd('create surface curve all')
fidesys.cmd('sweep surface 1 perpendicular distance 0.005 keep')
fidesys.cmd('delete surface 1')
fidesys.cmd('volume all size auto factor 3')
fidesys.cmd('mesh volume all')
fidesys.cmd('create material 1 ')
fidesys.cmd('modify material 1 name \'alum6061-T6\'')
fidesys.cmd('modify material 1 set property \'MODULUS\' value 6.894e+10')
fidesys.cmd('modify material 1 set property \'POISSON\' value 0.3')
fidesys.cmd('modify material 1 set property \'DENSITY\' value 2849')
fidesys.cmd('block 1 add volume all')
fidesys.cmd('block 1 cs 1 element solid order 1')
fidesys.cmd('block 1 material 1 cs 1 element solid order 1')
fidesys.cmd('create displacement  on surface 13  dof all fix 0 ')
fidesys.cmd('create displacement  on surface 7  dof 2 dof 3 fix 0 ')
fidesys.cmd('create distributed force on surface 7  force value 2500 moment value 0 direction 1 0 0 equivalent')
fidesys.cmd('analysis type spectrum elasticity dim3 preload on')
fidesys.cmd('spectrum type random_vibration')
fidesys.cmd('eigenvalue find 6 smallest')
fidesys.cmd(f"calculation start path '{BASE_DIR}/random_vibration.pvd'")
nodalCoords = [] # массив координат узлов
nodesId = [] # массив id узлов
elementsNodesConnect = [] # массив связи узлов и элементов
nodeCount = cubit.get_node_count() # количество узлов
elementCount = cubit.get_element_count() # количество элементов
for nodeId in range(1,nodeCount+1):
    nodesId.append(nodeId) # заполнение массива id узлов
    x,y,z = cubit.get_nodal_coordinates(nodeId) # получение координат для каждого узла
    nodalCoords.append([x,y,z]) # заполнение массива координат узлов
for elementId in range(1,elementCount+1):
     nodeIdList = cubit.get_expanded_connectivity("hex", elementId)  #поиск кортежа узлов для каждого гексаэдрального элемента
     nodeIdList = list(nodeIdList)   # конвертация кортежа в массив 
     elementsNodesConnect.append(list(map(lambda el: el-1, nodeIdList)))  # в каждом массиве узлов от каждого элемента необходимо отнять единицу
fc.delete_application() # удаление задачи из памяти

# блок для работы с Viewer
reader = OpenDataFile(f'{BASE_DIR}/random_vibration.pvd') # создание ридера
dataReader = servermanager.Fetch(reader) # получение всех данных с ридера
fieldData = dataReader.GetFieldData()  # получение всех полей из данных
freqField = fieldData.GetArray('Eigen Values') # получение массива частот
statisticalFreq = vtk_to_numpy(freqField)[0][0] # статистическая частота как первая собственная частота модели
pointData = dataReader.GetPointData() # получение всех значений точек из данных
arrayNodeId = pointData.GetArray('Node ID') # получение id узлов
nodeIdNumpy = vtk_to_numpy(arrayNodeId)  # получение numpy массива id узлов
randomVibrationReader = RandomVibration(reader) # получение ридера для фильтра Random Vibration
randomVibrationReader.Valuesforprocessing = ['Stress'] # задание данных для вывода
randomVibrationReader.Numberoffrequencysamples = 200 # получение количества точек для вывода спектральной плотности мощности массивов (можно оставить по дефолту, так как в этой задаче это не пригодится) 
randomVibrationReader.ConstantModeDampingRatio = 0.02 # задание коэффициентов демпфирования
randomVibrationReader.InputPSDdata = f'{BASE_DIR}/accel.csv' # задание кривой ускорения
randomVibrationReader.UpdatePipeline() # обновление данных для текущего шага
randomVibrationOutput_1 = randomVibrationReader.__getitem__(1) # второй вход для фильтра Random Vibration
randomVibrationOutput_1.UpdatePipeline() # обновление данных для текущего шага
extract = ExtractSelection(randomVibrationOutput_1) # Создание выборки для входа randomVibrationOutput_1
selection = _select('POINT',f'((NodeID <= {nodeCount}))',reader) # выбор всех точек модели
extract.Selection = selection # Заполние выборки всеми узлами
selectedData = servermanager.Fetch(extract) # получение всех данных с выборки
randomPointData = selectedData.GetPointData() # получение всех значений точек из выборки
arrayStressRMS = randomPointData.GetArray('Stress_RMS') # получение массива напряжений
stressRMSNumpy  = vtk_to_numpy(arrayStressRMS) # получение numpy массива напряжений
stressRMSNumpy_XX = stressRMSNumpy[:,0] # получение numpy массива напряжений XX
# Блок для подсчета массивов Life и Damage
stressArrayToWrite = [] # масив напряжений для записи в vtu
for id in range(1,nodeCount + 1):
    index = np.where(nodeIdNumpy==id)[0][0] # получение индекса узла в numpy массиве
    stressArrayToWrite.append(stressRMSNumpy_XX[index]) # заполнение массива напряжений для записи
Damage = [] # массив коэффициента разрушаемости
denominatorLife = [] # массив для знаменателя выражения долговечности
n1sigmaPersecond = 0.683*statisticalFreq # количество циклов  для размахов меньше 1sigma напряжений в единицу времени
n2sigmaPersecond = 0.271*statisticalFreq # количество циклов  для размахов меньше 2sigma напряжений в единицу времени
n3sigmaPersecond = 0.0433 * statisticalFreq # количество циклов  для размахов меньше 3sigma напряжений в единицу времени
N1sigma = list(map(lambda el: 10**((math.log10(el)- B)/K), stressArrayToWrite))   # количество допускаемых циклов до разрушения  для размахов меньше 1sigma напряжений 
N2sigma = list(map(lambda el: 10**((math.log10(2*el)- B)/K), stressArrayToWrite)) # количество допускаемых циклов до разрушения  для размахов меньше 2sigma напряжений 
N3sigma = list(map(lambda el: 10**((math.log10(3*el)- B)/K), stressArrayToWrite)) # количество допускаемых циклов до разрушения  для размахов меньше 3sigma напряжений 
n1sigma =n1sigmaPersecond * DURATION_TIME   # количество циклов  для размахов меньше 1sigma в течение всего испытания
n2sigma = n2sigmaPersecond * DURATION_TIME  # количество циклов  для размахов меньше 1sigma в течение всего испытания
n3sigma = n3sigmaPersecond * DURATION_TIME  # количество циклов  для размахов меньше 1sigma в течение всего испытания
Damage1sigma =  list(map(lambda el: n1sigma/el, N1sigma))  # коэффициент разрушаемости для   для размахов меньше 1sigma
Damage2sigma =  list(map(lambda el: n2sigma/el, N2sigma))   # коэффициент разрушаемости для   для размахов меньше 2sigma
Damage3sigma =  list(map(lambda el: n3sigma/el, N3sigma))   # коэффициент разрушаемости для   для размахов меньше 3sigma
denominatorLife1sigma = list(map(lambda el: n1sigmaPersecond/el, N1sigma))  # массив для знаменателя выражения долговечности  1sigma
denominatorLife2sigma = list(map(lambda el: n2sigmaPersecond/el, N2sigma))  # массив для знаменателя выражения долговечности  2sigma
denominatorLife3sigma = list(map(lambda el: n3sigmaPersecond/el, N3sigma))  # массив для знаменателя выражения долговечности  3sigma
for i in range(len(stressArrayToWrite)):
    Damage.append(Damage1sigma[i] +  Damage2sigma[i] + Damage3sigma[i]) # заполнение массива коэффициента повреждаемости
    denominatorLife.append(denominatorLife1sigma[i] + denominatorLife2sigma[i] + denominatorLife3sigma[i]) # заполнение массива  знаменателя выражения долговечности
Life = list(map(lambda el: el**-1, denominatorLife))  # получение массива долговечности (в секундах)
Life = list(map(lambda el: el/DURATION_TIME/365, Life)) # конвертация массива долговечности из секунд в года
# блок для записи vtu - файла
cells = [
    ("hexahedron", elementsNodesConnect) # ячейки в виде тетраэдров первого порядка
]
mesh = meshio.Mesh(
    nodalCoords,
    cells,
    point_data= {
        "NodeId":nodesId,
        "Stress_RMS":stressArrayToWrite,
        "Damage":Damage,
        "Life":Life
        }
) # формирование сетки и значений в узлах заданных массивов
mesh.write(
    os.path.join(BASE_DIR ,"fatique.vtu")  
) # запись vtu - файла

Запустите командную строку из папки. Затем создайте виртуальную среду (например, venv). Пропишите в командной строке:

python -m virtualenv venv

Запустите виртуальную среду командой:

.\venv\Scripts\activate

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

Установите в окружающую среду библиотеки meshio, numpy (версия < 1.24):

pip install meshio
pip install numpy==1.23.5

Выйти из виртуальной среды:

deactivate

Запустите файл calculate.py. Результатом должен быть собранный vtu – файл “fatique”, который успешно открывается с помощью FidesysViewer. В данном файле должны присутствовать следующие значения точек (node_id, Stress_RMS, Damage, Life)

Анализ результатов

Откройте файл с результатами. Появится окно FidesysViewer, в котором вы сможете ознакомиться с результатами расчёта.

На верхней панели выберите данные результата расчета для отображения - Damage. Для корректности изображения эпюры необходимо перевести эпюры в логарифмический масштаб.

Нажмите Визуализация и закройте окно. После чего растяните высоту шкалы.

Затем на верхней панели выберите данные результаты расчета для отображения - Life. Для корректности изображения эпюры необходимо перевести эпюры в логарифмический масштаб.

Так как для поля “Срок службы” опасное значение является минимальным, то его необходимо покрасить в красный цвет. Для этого в “редактировать свойства цифровой легенды” необходимо нажать “Инвертировать передаточные соотношения”.

Нажмите Визуализация и закройте окно. После чего растяните высоту шкалы.

Таким образом, при воздействии на образец широкополостной вибрации в течение 24 часа коэффициент повреждаемости меньше 1, следовательно данный образец выдержит испытание на усталость. При этом минимальный срок службы образца составит 25 лет.