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

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

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

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

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

importGeometry(model,'Gen3Shoulder.stl');

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

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

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

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  

Лучший способ визуализировать формы режима - это анимировать гармоническое движение на соответствующих частотах. The 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.

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 смотрите раздел Pressure Function в нижней части этой страницы. Эта функция не имеет явной зависимости от частоты. Поэтому, в частотный диапазон, эта нагрузка давления действует на все частоты решения.

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