В этом примере показано, как включать затухание в анализ переходных процессов простого консольного луча.
Модель затухания является основным вязким затуханием, распределенным однородно через объем луча. Луч деформирован путем применения внешней загрузки в совете луча и затем выпущен во время . Этот пример не использует дополнительной загрузки, таким образом, смещение уменьшений луча как функция времени из-за затухания. Пример использует плоское напряжение модальные, статические модели, и анализа переходных процессов в его рабочем процессе с тремя шагами:
Выполните модальный анализ, чтобы вычислить основную частоту луча и ускорить расчеты для анализа переходных процессов.
Найдите, что статическое решение луча с вертикальной загрузкой в совете использует в качестве начального условия для переходной модели.
Выполните анализ переходных процессов с и без затухания.
Затухание обычно выражается как процент критического затухания, , для выбранной частоты вибрации. Этот пример использует , который составляет три процента критического затухания.
Пример задает значения параметров с помощью имперской системы модулей. Можно заменить их на значения, заданные в метрической системе. Если вы делаете так, гарантируете, что задаете все значения в примере с помощью той же системы.
Создайте модальную аналитическую модель для проблемы плоского напряжения.
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: [8×1 double] ModeShapes: [1×1 FEStruct] Mesh: [1×1 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
Луч деформирован путем применения внешней загрузки в ее совете и затем выпущен во время . Найдите начальное условие для анализа переходных процессов при помощи статического решения луча с вертикальной загрузкой в совете.
Создайте статическую модель плоского напряжения.
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: [6511×2 double] InitialVelocity: [6511×2 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-смещение в совете является синусоидальной функцией времени с амплитудой экспоненциально descreasing со временем.
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')