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

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

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'

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

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=0$$. Найдите начальное условие для анализа переходных процессов при помощи статического решения луча с вертикальной загрузкой в совете.

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

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

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

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');

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

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')