Модальный анализ камертона

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

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

model = createpde('structural','modal-solid');

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

importGeometry(model,'TuningFork.stl');
figure
pdegplot(model)

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

E = 210E9;
nu = 0.3;
rho = 8000;
structuralProperties(model,'YoungsModulus',E, ...
                           'PoissonsRatio',nu, ...
                           'MassDensity',rho);

Сгенерируйте mesh.

generateMesh(model,'Hmax',0.001);

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

RF = solve(model,'FrequencyRange',[-1,4000]*2*pi);

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

modeID = 1:numel(RF.NaturalFrequencies);

Выразите получившиеся частоты в Гц путем деления их на $2\pi$. Отображение частот в таблице.

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

    Mode    Frequency
    ____    _________

      1     0.0072398
      2     0.0033543
      3     0.0025636
      4     0.0039618
      5     0.0053295
      6     0.0094544
      7        460.42
      8        706.34
      9        1911.5
     10        2105.5
     11        2906.5
     12        3814.7

Поскольку в этом примере нет граничных ограничений, модальные результаты включают твердые режимы тела. Первые шесть частот, близких к нулю, указывают на шесть режимов твёрдого тела 3-D твердого тела. Первый гибкий режим является седьмым режимом с частотой около 460 Гц.

Лучший способ визуализировать формы режима - это анимировать гармоническое движение на их соответствующих частотах. The animateSixTuningForkModes функция анимирует шесть гибких режимов, которые являются режимами с 7 по 12 в модальных результатах RF.

frames  = animateSixTuningForkModes(RF);

Для воспроизведения анимации используйте следующую команду:

movie(figure('units','normalized','outerposition',[0 0 1 1]),frames,5,30)

В первом режиме две осциллирующие стойки камертона балансируют поперечные силы на указателе. Следующий режим с этим эффектом - пятый гибкий режим с частотой 2906,5 Гц. Эта частота примерно в 6,25 раза больше основной частоты 460 Гц.

Переходный анализ камертона

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

Во-первых, создайте модель структурного переходного анализа.

tmodel = createpde('structural','transient-solid');

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

importGeometry(tmodel,'TuningFork.stl');

Сгенерируйте mesh.

mesh = generateMesh(tmodel,'Hmax',0.005);

Задайте модуль Юнга, отношение Пуассона и массовую плотность.

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

Идентифицируйте грани для применения граничных ограничений и нагрузок путем построения графика геометрии с метками граней.

figure('units','normalized','outerposition',[0 0 1 1])
pdegplot(tmodel,'FaceLabels','on')
view(-50,15)
title 'Geometry with Face Labels'

Наложите достаточные граничные ограничения, чтобы предотвратить движение твердого тела при приложенной загрузке. Обычно вы держите камертон вручную или монтируете его на таблице. Упрощенным приближением к этому граничному условию является фиксация области вблизи пересечения выступов и указателя (граней 21 и 22).

structuralBC(tmodel,'Face',[21,22],'Constraint','fixed');

Аппроксимируйте импульсную загрузку на поверхность стойки, прикладывая нагрузку под давлением в течение очень малой части периода времени основного режима. При помощи этого очень короткого импульса давления вы гарантируете, что только основной режим камертона возбужден. Чтобы вычислить период времени T основного режима, использовать результаты модального анализа.

T = 2*pi/RF.NaturalFrequencies(7);

Задайте загрузку давлением на стойку как короткий прямоугольный импульс давления.

structuralBoundaryLoad(tmodel,'Face',11,'Pressure',5E6,'EndTime',T/300);

Примените нулевое перемещение и скорость в качестве начальных условий.

structuralIC(tmodel,'Displacement',[0;0;0],'Velocity',[0;0;0]);

Решите переходную модель на 50 периодов основного режима. Дискретизируйте динамику 60 раз за период основного режима.

ncycle = 50;
samplingFrequency = 60/T;
tlist = linspace(0,ncycle*T,ncycle*T*samplingFrequency);
R = solve(tmodel,tlist)

R = 
  TransientStructuralResults with properties:

     Displacement: [1×1 FEStruct]
         Velocity: [1×1 FEStruct]
     Acceleration: [1×1 FEStruct]
    SolutionTimes: [1×3000 double]
             Mesh: [1×1 FEMesh]

Постройте график timeseries колебаний совета оголовка, который является поверхностью 12. Найдите узлы на грани совета и постройте график y-компонент перемещения с течением времени, с использованием одного из этих узлов.

excitedTineTipNodes = findNodes(mesh,'region','Face',12);
tipDisp = R.Displacement.uy(excitedTineTipNodes(1),:);

figure
plot(R.SolutionTimes,tipDisp)
title('Transverse Displacement at Tine Tip')
xlim([0,0.1])
xlabel('Time')
ylabel('Y-Displacement')

Выполните быстрое преобразование Фурье (FFT) на timeseries перемещения совета, чтобы увидеть, что частота вибрации камертона близка к его основной частоте. Небольшое отклонение от основной частоты, вычисленной в неограниченном модальном анализе, появляется из-за ограничений, наложенных в переходном анализе.

[fTip,PTip] = tuningForkFFT(tipDisp,samplingFrequency);
figure
plot(fTip,PTip) 
title({'Single-sided Amplitude Spectrum', 'of Tip Vibration'})
xlabel('f (Hz)')
ylabel('|P1(f)|')
xlim([0,4000])

Поперечная вибрация зубцов заставляет указатель колебаться в осевом направлении с той же частотой. Чтобы наблюдать эту вибрацию, постройте график временного timeseries осевого перемещения торцевой поверхности указателя.

baseNodes = tmodel.Mesh.findNodes('region','Face',6);
baseDisp = R.Displacement.ux(baseNodes(1),:);
figure
plot(R.SolutionTimes,baseDisp)
title('Axial Displacement at the End of Handle')
xlim([0,0.1])
ylabel('X-Displacement')
xlabel('Time')

Выполните БПФ timeseries осевой вибрации указателя. Эта частота вибрации также близка к своей основной частоте.

[fBase,PBase] = tuningForkFFT(baseDisp,samplingFrequency);
figure
plot(fBase,PBase) 
title({'Single-sided Amplitude Spectrum', 'of Base Vibration'})
xlabel('f (Hz)')
ylabel('|P1(f)|')
xlim([0,4000])