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

Предположим, что один конец роботизированной руки зафиксирован, найдите собственные частоты и формы режимов.

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

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

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

importGeometry(model,'Gen3Shoulder.stl');

Создайте сетку.

generateMesh(model);
pdemesh(model)

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

E = 1.5E11;          
nu = 0.3;           
rho = 2000;            
structuralProperties(model,'YoungsModulus',E, ...
                           'PoissonsRatio',nu, ...
                           'MassDensity',rho);

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

figure
pdegplot(model,'FaceLabels','on')
view([-1 2])
title('Shoulder Link Geometry with Face Labels')

Плечевое звено закреплено на одном конце (грань 3) и соединено с подвижным звеном на другом конце (грань 4). Применить фиксированное граничное условие к грани 3.

structuralBC(model,'Face',3,'Constraint','fixed');

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

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

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

modeID = 1:numel(RF.NaturalFrequencies);

Выразить результирующие частоты в Гц, разделив их на $\mathrm{}\mathrm{2\delta }$. Отображение частот в таблице.

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

    Mode    Frequency
    ____    _________

     1       1947.2  
     2         2662  
     3       4982.3  
     4       5112.6  
     5       7819.5  
     6       8037.1  
     7         9361  

Лучший способ визуализации форм моды - анимация гармонического движения на соответствующих частотах. animateSixLinkModes функция анимирует первые шесть режимов. Полученный график показывает области доминирующей деформации под нагрузкой.

figure
frames  = animateSixLinkModes(RF);

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

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

Анализ частотной характеристики

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

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

fmodel = createpde('structural','frequency-solid');

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

importGeometry(fmodel,'Gen3Shoulder.stl');

Создайте сетку.

mesh = generateMesh(fmodel);

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

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

Плечевое звено закреплено на одном конце (грань 3) и соединено с подвижным звеном на другом конце (грань 4). Применить фиксированное граничное условие к грани 3.

structuralBC(fmodel,'Face',3,'Constraint','fixed');

Оцените давление, которое подвижное звено прикладывает к грани 4, когда рычаг несет нагрузку. На этом рисунке показаны две половины грани 4, разделенные в центре по координате Y.

Используйте pressFcnFR для приложения граничной нагрузки к грани 4. Эта функция применяет сигналы давления нажатия и скручивания. Компонент давления толкания является однородным. Компонент скручивания оказывает положительное давление на левую сторону и отрицательное давление на правую сторону грани. Для определения pressFcnFR см. раздел «Функция давления» в нижней части этой страницы. Эта функция не имеет явной зависимости от частоты. Поэтому в частотной области эта нагрузка давления действует на все частоты решения.

structuralBoundaryLoad(fmodel,'Face',4,'Pressure',@(region,state) pressFcnFR(region, state),'Vectorized','on');

Определите список частот для решения от 0 до 3500 Гц с 200 шагами.

flist = linspace(0,3500,200)*2*pi;

Решите модель с помощью модального решателя частотных характеристик, указав объект RF результатов модели в качестве одного из входов.

R = solve(fmodel,flist,'ModalResults',RF)

R = 
  FrequencyStructuralResults with properties:

           Displacement: [1x1 FEStruct]
               Velocity: [1x1 FEStruct]
           Acceleration: [1x1 FEStruct]
    SolutionFrequencies: [1x200 double]
                   Mesh: [1x1 FEMesh]

Постройте график частотной характеристики в точке на загруженной грани. Точка на грани 4, расположенная при максимальной нагрузке отрицательного давления (0.003; 0.0436; 0.1307). Выполните интерполяцию смещения до этой точки и постройте график результата.

queryPoint  = [0.003; 0.0436; 0.1307];
queryPointDisp = R.interpolateDisplacement(queryPoint); 

figure
plot(R.SolutionFrequencies/2/pi,abs(queryPointDisp.uy))
title('Transverse Displacement at Point on Loaded Face')
xlabel('Frequency (Hz)')
ylabel('Y-Displacement')
xlim([0.0000 3500])

Пик отклика происходит вблизи 2662 Гц, что близко ко второму режиму вибрации. Меньший отклик также возникает при первом режиме, близком к 1947 Гц.

Найдите индекс пиковой частоты отклика с помощью max с двумя выходными аргументами. Второй выходной аргумент обеспечивает индекс пиковой частоты.

[M, I] = max(abs(queryPointDisp.uy))

M = 1.1256e-04

I = 152

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

RD = struct();
RD.ux = R.Displacement.ux(:,I);
RD.uy = R.Displacement.uy(:,I);
RD.uz = R.Displacement.uz(:,I);

figure('units','normalized','outerposition',[0 0 1 1]);

subplot(2,2,1)
pdeplot3D(fmodel,'ColorMapData',R.Displacement.ux(:,I), ...
                 'Deformation',RD,'DeformationScaleFactor',1);
title('X-Displacement')

subplot(2,2,2)
pdeplot3D(fmodel,'ColorMapData',R.Displacement.uy(:,I), ...
                 'Deformation',RD,'DeformationScaleFactor',1);
title('Y-Displacement')

subplot(2,2,3)
pdeplot3D(fmodel,'ColorMapData',R.Displacement.uz(:,I), ...
                 'Deformation',RD,'DeformationScaleFactor',1);
title('Z-Displacement')

subplot(2,2,4)
pdeplot3D(fmodel,'ColorMapData',R.Displacement.Magnitude(:,I), ...
                 'Deformation',RD,'DeformationScaleFactor',1);
title('Magnitude')

Функция давления

Определите функцию давления, pressFcnFR, чтобы вычислить сигналы давления толчка и скручивания. Компонент давления толкания является однородным. Компонент давления скручивания оказывает положительное давление на левую сторону и отрицательное давление на правую сторону грани. Значение нагрузки давления скручивания увеличивается в параболическом распределении от минимального в точке C к положительному пику при L и к отрицательному пику при R. Коэффициент давления скручивания для параболического распределения, полученного в pressFcnFR умножается на синусоидальную функцию с величиной 0,1 МПа. Равномерное давление толкания составляет 10 кПа.

function p = pressFcnFR(region,~)

meanY = mean(region.y);
absMaxY = max(abs(region.y));
scaleFactor = zeros(size(region.y));

% Find IDs of the points on the left and right halves of the face
% using y-coordinate values.
leftHalfIdx = region.y <= meanY;
rightHalfIdx = region.y >= meanY;

% Define a parabolic scale factor for each half of the face.
scaleFactor(leftHalfIdx) = ((region.y(leftHalfIdx) - meanY)/absMaxY).^2;
scaleFactor(rightHalfIdx) = -((region.y(rightHalfIdx) - meanY)/absMaxY).^2;

p = 10E3 + 0.1E6*scaleFactor;

end