Вибрация квадратного диска

Этот пример показывает, как вычислить режимы вибрации и частоты 3-D просто поддерживаемого квадратного упругого диска.

Свойства размерностей и материала пластины взяты из стандартной задачи бенчмарка конечных элементов, опубликованной NAFEMS, FV52 (См. Ссылку).

Во-первых, создайте контейнер структурной модели для вашей 3-D задачи модального анализа. Это контейнер, который содержит геометрию, свойства материала, нагрузки на тело, краевые нагрузки, граничные ограничения и mesh.

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

Импорт STL- файла модели простого диска с помощью importGeometry функция. Эта функция восстанавливает грани, ребра и вершины модели. Он может объединять некоторые грани и ребра, поэтому числа могут отличаться от чисел родительской модели CAD.

importGeometry(model,'Plate10x10x1.stl');

Постройте график геометрии и включите метки граней. Метки граней нужны при определении граничных условий.

figure
hc = pdegplot(model,'FaceLabels','on');
hc(1).FaceAlpha = 0.5;
title('Plate with Face Labels')

Задайте модуль упругости стали, коэффициент Пуассона и плотность материала.

structuralProperties(model,'YoungsModulus',200e9, ...
                           'PoissonsRatio',0.3, ...
                           'MassDensity',8000);

В этом примере единственным граничным условием является нуль $$z$$-размещение на четырех гранях ребер. Эти ребра имеют метки с 1 по 4.

structuralBC(model,'Face',1:4,'ZDisplacement',0);

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

generateMesh(model,'Hmin',1.3);
figure 
pdeplot3D(model);
title('Mesh with Quadratic Tetrahedral Elements');

Для сравнения с опубликованными значениями загрузите опорные частоты в Гц.

refFreqHz = [0 0 0 45.897 109.44 109.44 167.89 193.59 206.19 206.19];

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

maxFreq = 1.1*refFreqHz(end)*2*pi;
result = solve(model,'FrequencyRange',[-0.1 maxFreq]);

Вычислите частоты в Гц.

freqHz = result.NaturalFrequencies/(2*pi);

Сравните базовую и вычисленную частоты (в Гц) для самых низких 10 режимов. Самые низкие трехмодовые формы соответствуют движению жесткого тела пластины. Их частоты близки к нулю.

tfreqHz = table(refFreqHz.',freqHz(1:10));
tfreqHz.Properties.VariableNames = {'Reference','Computed'};
disp(tfreqHz);

    Reference     Computed 
    _________    __________

          0      2.8632e-05
          0      7.4852e-06
          0      2.0318e-05
     45.897          44.871
     109.44          109.74
     109.44          109.77
     167.89          168.59
     193.59          193.74
     206.19          207.51
     206.19          207.52

Вы видите хорошее согласие между вычисляемой и опубликованной частотами.

Постройте график третьего компонента ($$z$$-компонент) решения для семи самых низких ненулевых частотных режимов.

h = figure;
h.Position = [100,100,900,600];
numToPrint = min(length(freqHz),length(refFreqHz));
for i = 4:numToPrint
    subplot(4,2,i-3);
    pdeplot3D(model,'ColorMapData',result.ModeShapes.uz(:,i));
    axis equal
    title(sprintf(['Mode=%d, z-displacement\n', ...
    'Frequency(Hz): Ref=%g FEM=%g'], ...
    i,refFreqHz(i),freqHz(i)));
end