Модальный анализ

Создайте модель модального анализа для задачи плоское напряжение.

modelM = createpde('structural','modal-planestress');

Создайте геометрию и включите ее в модель. Предположим, балка имеет длину 5 дюймов и толщину 0,1 дюйма.

width = 5; 
height = 0.1;

gdm = [3;4;0;width;width;0;0;0;height;height];
g = decsg(gdm,'S1',('S1')');
geometryFromEdges(modelM,g);

Постройте график геометрии с метками кромок.

figure; 
pdegplot(modelM,'EdgeLabels','on'); 
axis equal
title 'Geometry With Edge Labels Displayed'

Определите максимальный размер элемента таким образом, чтобы толщина балки составляла пять элементов. Создайте сетку.

hmax = height/5;
msh = generateMesh(modelM,'Hmax',hmax);

Укажите модуль Юнга, коэффициент Пуассона и массовую плотность стали.

E = 3.0e7; 
nu = 0.3; 
rho = 0.3/386;
structuralProperties(modelM,'YoungsModulus',E, ...
                            'PoissonsRatio',nu, ...
                            'MassDensity',rho);

Укажите, что левая кромка балки является фиксированной границей.

structuralBC(modelM,'Edge',4,'Constraint','fixed');

Решить проблему для диапазона частот от 0 кому 1e5. Рекомендуется использовать значение, которое немного меньше ожидаемой самой низкой частоты. Таким образом, использование -0.1 вместо 0.

res = solve(modelM,'FrequencyRange',[-0.1,1e5]')

res = 
  ModalStructuralResults with properties:

    NaturalFrequencies: [8x1 double]
            ModeShapes: [1x1 FEStruct]
                  Mesh: [1x1 FEMesh]

По умолчанию решатель возвращает круговые частоты.

modeID = 1:numel(res.NaturalFrequencies);

Выразить результирующие частоты в Гц, разделив их на 2π. Отображение частот в таблице.

tmodalResults = table(modeID.',res.NaturalFrequencies/(2*pi));
tmodalResults.Properties.VariableNames = {'Mode','Frequency'};
disp(tmodalResults)

    Mode    Frequency
    ____    _________

     1       126.94  
     2       794.05  
     3       2216.8  
     4       4325.3  
     5       7110.7  
     6       9825.9  
     7        10551  
     8        14623  

Вычислите аналитическую основную частоту (Гц), используя теорию луча.

I = height^3/12;
freqAnalytical = 3.516*sqrt(E*I/(width^4*rho*height))/(2*pi)

freqAnalytical = 126.9498

Сравните аналитический результат с числовым результатом.

freqNumerical = res.NaturalFrequencies(1)/(2*pi)

freqNumerical = 126.9416

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

longestPeriod = 1/freqNumerical

longestPeriod = 0.0079

Постройте график y-составляющей решения для самой низкой частоты луча.

figure;
pdeplot(modelM,'XYData',res.ModeShapes.uy(:,1))
title('Lowest Frequency Vibration Mode')
axis equal

Начальное смещение из статического решения

Балка деформируется путем приложения внешней нагрузки на ее кончике и затем освобождается в момент времени $$t\mathrm{=}$$0. Найдите исходное условие для анализа переходных процессов, используя статическое решение балки с вертикальной нагрузкой на наконечнике.

Создание статической модели «плоскость-напряжение».

modelS = createpde('structural','static-planestress');

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

geometryFromEdges(modelS,g);
modelS.Mesh = msh;

Задайте одинаковые значения для модуля Юнга, коэффициента Пуассона и плотности массы материала.

structuralProperties(modelS,'YoungsModulus',E, ...
                           'PoissonsRatio',nu, ...
                           'MassDensity',rho);
                       

Задайте такую же зависимость на левом конце балки.

structuralBC(modelS,'Edge',4,'Constraint','fixed');

Примените статическую вертикальную нагрузку с правой стороны балки.

structuralBoundaryLoad(modelS,'Edge',2,'SurfaceTraction',[0;1]);

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

Rstatic = solve(modelS);

Анализ переходных процессов

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

Создайте модель «переходная плоскость-напряжение».

modelT = createpde('structural','transient-planestress');

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

geometryFromEdges(modelT,g);
modelT.Mesh = msh;

Задайте одинаковые значения для модуля Юнга, коэффициента Пуассона и плотности массы материала.

structuralProperties(modelT,'YoungsModulus',E, ...
                            'PoissonsRatio',nu, ...
                            'MassDensity',rho);
                       

Задайте такую же зависимость на левом конце балки.

structuralBC(modelT,'Edge',4,'Constraint','fixed');

Укажите начальное условие с помощью статического решения.

structuralIC(modelT,Rstatic)

ans = 
  NodalStructuralICs with properties:

    InitialDisplacement: [6511x2 double]
        InitialVelocity: [6511x2 double]

Решите неразвернутую модель переходного процесса для трех полных периодов, соответствующих самому низкому режиму вибрации.

tlist = 0:longestPeriod/100:3*longestPeriod;
resT = solve(modelT,tlist,'ModalResults',res);

Интерполировать смещение в конце балки.

intrpUt = interpolateDisplacement(resT,[5;0.05]);

Смещение на кончике является синусоидальной функцией времени с амплитудой, равной начальному y-смещению. Этот результат согласуется с решением простой пружинно-массовой системы.

plot(resT.SolutionTimes,intrpUt.uy)
grid on
title('Undamped Solution')
xlabel('Time')
ylabel('Tip of beam displacement')

Теперь решите модель с демпфированием, равным 3% критического демпфирования.

zeta = 0.03;
omega = 2*pi*freqNumerical;
structuralDamping(modelT,'Zeta',zeta);
resT = solve(modelT,tlist,'ModalResults',res);

Интерполировать смещение в конце балки.

intrpUt = interpolateDisplacement(resT,[5;0.05]);

Y-смещение на кончике является синусоидальной функцией времени с экспоненциально уменьшающейся со временем амплитудой.

figure
hold on
plot(resT.SolutionTimes,intrpUt.uy)
plot(tlist,intrpUt.uy(1)*exp(-zeta*omega*tlist),'Color','r')
grid on
title('Damped Solution')
xlabel('Time')
ylabel('Tip of beam displacement')