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

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

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