Динамика ослабленного консольного луча

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

Луч моделируется с плоской формулировкой эластичности напряжения. Модель затухания является основным вязким затуханием, распределенным однородно через объем луча. Несколько анализов переходных процессов выполняются, где луч деформирован в начальную форму и затем выпущен во время, $t = 0$ . Исследования с и без затухания рассматриваются. Две начальных формы смещения рассматриваются. В первом луч деформирован в форму, соответствующую самому низкому режиму вибрации. Во втором луч деформирован путем применения внешней загрузки в совете луча. Никакая дополнительная загрузка не применяется в этом примере так, в ослабленных случаях, смещении затуханий луча как функция времени из-за затухания.

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

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

$D = α M + β K$

Распространено выразить затухание как процент критического затухания, $ξ$ , для выбранной частоты вибрации. Для данной частоты, $ω_{i}$ , следующее выражение имеет отношение $ξ$ к $α$ и $β$ .

$ξ = \frac{α}{2 ω_{i}} + \frac{β ω_{i}}{2}$

В этом примере мы зададим $ξ = 0 . 03$ (три процента критического затухания) и $β$ равняйтесь нулю так, чтобы $α$ может быть вычислен как $α = 2 \times . 03 ω_{i}$

Этот пример задает значения параметров с помощью имперской системы модулей. Можно заменить их на значения, заданные в метрической системе. Если вы делаете так, гарантируете, что задаете все значения в примере с помощью той же системы.

Излучите параметры

Задайте размерности луча. Это 5 дюймов длиной и 0,1 дюйма толщиной.

width = 5; 
height = 0.1;

Укажите, что материал является сталью.

E = 3.0e7; 
nu = 0.3; 
rho = 0.3/386;

Вычислите матрицу коэффициентов от свойств материала.

G = E/(2.*(1+nu));
mu = 2.0*G*nu/(1-nu);

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

I = height^3/12; 
A = height;

eigValAnalytical = 1.8751^4*E*I/(rho*A*width^4); 
freqAnalytical = sqrt(eigValAnalytical)/(2*pi);

Setup задач

Создайте PDEModel с двумя независимыми переменными, чтобы представлять анализ.

numberOfPDE = 2;
model = createpde(numberOfPDE);

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

gdm = [3;4;0;width;width;0;0;0;height;height];
g = decsg(gdm,'S1',('S1')');

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

figure; 
pdegplot(g,'EdgeLabels','on'); 
axis equal
title 'Geometry With Edge Labels Displayed'

Предоставьте модели определение геометрии.

geometryFromEdges(model,g);

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

c = [2*G+mu;0;G;0;G;mu;0;G;0;2*G+mu];
f = [0;0];
a = 0;
m = rho;
specifyCoefficients(model,'m',m,'d',0,'c',c,'a',a,'f',f);

Задайте следующее граничное условие, чтобы зафиксировать (смещения равный нуль) левое ребро луча.

applyBoundaryCondition(model,'dirichlet','Edge',4,'u',[0 0]);

Задайте максимальный размер элемента (5 элементов через толщину луча) и сгенерируйте mesh.

hmax = height/5;
msh=generateMesh(model,'Hmax',hmax,'MesherVersion','R2013a');

Вычисление режимов вибрации и частот

Используйте solvepdeeig и затем вычислите режим вибрации самой низкой частоты.

res = solvepdeeig(model, [0,1e6]');

              Basis= 10,  Time=   0.47,  New conv eig=  3
End of sweep: Basis= 10,  Time=   0.47,  New conv eig=  3
              Basis= 13,  Time=   0.83,  New conv eig=  1
End of sweep: Basis= 13,  Time=   0.83,  New conv eig=  1

eigenVec = res.Eigenvectors;
eigenVal = res.Eigenvalues;

Цветной график y- смещение режима вибрации самой низкой частоты.

freqNumerical = sqrt(eigenVal(1))./(2*pi);
longestPeriod = 1/freqNumerical;

Постройте деформированную форму луча со смещениями, масштабируемыми случайным фактором.

scaleFactor = 20;
figure;
[p,e,t] = meshToPet(msh);
pdeplot(p+scaleFactor*eigenVec',e,t,'XYData',real(eigenVec(:,2)));
title('Lowest Frequency Vibration Mode');
axis equal
xlabel('Inches');
ylabel('Inches');

drawnow

fprintf('Lowest beam frequency (Hz). Analytical = %8.3e, Numerical = %8.3e\n', ...
  freqAnalytical,freqNumerical);

Lowest beam frequency (Hz). Analytical = 1.269e+02, Numerical = 1.269e+02

Анализ переходных процессов, начальное смещение от первой формы режима

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

Для удобства мы будем масштабировать форму собственного вектора так, чтобы y- смещение в совете составляет.1 дюйм. Это делает переходные графики истории более простыми.

uEigenTip = eigenVec(2,2);
u0TipDisp = .1;

u0 = u0TipDisp/uEigenTip*eigenVec;

Сначала решите незатухающую систему.

Вычислите решение в течение трех полных периодов.

tlist = 0:longestPeriod/100:3*longestPeriod;

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

R = createPDEResults(model, u0(:));
ice = icEvaluator(R);

Установите начальные условия и решите систему.

setInitialConditions(model, @ice.computeIC, 0);
tres = solvepde(model,tlist);

Смещение в совете является синусоидальной функцией времени с амплитудой, равной начальному y-смещению. Это соглашается с решением простой пружинно-массовой системы.

titl = 'Initial Displacements from Lowest Eigenvector, Undamped Solution';
cantileverBeamTransientPlot(tres,titl);

Теперь решите систему с затуханием равного 3% критических для этой самой низкой частоты вибрации.

fem = assembleFEMatrices(model);
zeta = .03; 
omega = 2*pi*freqNumerical; 
alpha = 2*zeta*omega;
dampmat = alpha*fem.M;
specifyCoefficients(model,'m',m,'d',dampmat,'c',c,'a',a,'f',f);
tres = solvepde(model,tlist);

Этот рисунок показывает y- смещение в совете как функция времени. Наложенный на этот график вторая кривая, которая показывает конверт максимальной амплитуды как функция времени, вычисленного от решения до одного ОДУ степени свободы. Как ожидалось решение для УЧП соглашается хорошо с аналитическим решением.

titl = 'Initial Displacements from Lowest Eigenvector, Damped Solution';
cantileverBeamTransientPlot(tres,titl);
hold on;
plot(tlist,u0TipDisp*exp(-zeta*omega*tlist),'Color','r');
legend('PDE','ODE Amplitude Decay','Location','southeast');

Анализ переходных процессов, начальное смещение из статического решения

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

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

P = 1.0;
pdeTipLoad = createpde(2);
pg = geometryFromEdges(pdeTipLoad,g);

Задайте уравнение, чтобы решить.

specifyCoefficients(pdeTipLoad,'m',0,'d',0,'c',c,'a',a,'f',f);
tipLoad = applyBoundaryCondition(pdeTipLoad,'neumann','Edge',2,'g',[0 P/height]);
clampedEdge = applyBoundaryCondition(pdeTipLoad,'dirichlet','Edge',4,'u',[0,0]);
msh=generateMesh(pdeTipLoad,'Hmax',hmax,'MesherVersion','R2013a');

statres = solvepde(pdeTipLoad);

Чтобы сделать сравнение со случаем собственного вектора более ясным, мы будем также масштабировать это статическое решение так, чтобы максимальный y- смещение в совете равняется.1 дюйму.

u = statres.NodalSolution;
uEigenTip = u(2,2); 
u0TipDisp = 0.1;
u0 = u0TipDisp/uEigenTip*u;

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

specifyCoefficients(model, 'm', m, 'd', 0, 'c', c, 'a', a, 'f', f);

Установите начальные условия и решите систему.

R = createPDEResults(model, u0(:));
ice = icEvaluator(R);
setInitialConditions(model, @ice.computeIC, 0);
tres = solvepde(model,tlist);

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

titl = 'Initial Displacements from Static Solution, Undamped Solution';
cantileverBeamTransientPlot(tres,titl);

Вычислите ослабленное решение с начальным смещением от статического анализа. Матрица затухания эквивалентна, который использовал в случае собственного вектора, выше.

specifyCoefficients(model, 'm', m, 'd', dampmat, 'c', c, 'a', a, 'f', f);
tres = solvepde(model,tlist);

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

titl = 'Initial Displacements from Static Solution, Damped Solution';
cantileverBeamTransientPlot(tres,titl);
hold on;
plot(tlist,u0TipDisp*exp(-zeta*omega*tlist),'Color','r');
legend('PDE','ODE Amplitude Decay','Location','southeast');

Служебная функция построения графика

Служебная функция для создания графиков y-смещения совета как функция времени.

type cantileverBeamTransientPlot.m

function cantileverBeamTransientPlot( tdres, titl )
%CANTILEVERBEAMTRANSIENTPLOT Plot y-displacement at the beam tip
%   tdres - Time-dependent results object representing displacements as a function of time
%   titl plot title


% Copyright 2013-2015 The MathWorks, Inc.

tlist = tdres.SolutionTimes;
uu = tdres.NodalSolution;
utip = uu(2,2,:);
figure; plot(tlist, utip(:)); grid on;
title(titl); xlabel('Time, seconds');
ylabel('Y-displacement at beam tip, inches');
drawnow;
end

Указатель на функцию для определения Начального условия (IC)

Функция оценки возвращает значение IC в любой точке в mesh. Объект результатов представляет значения и функцию интерполяции, которую он обеспечил, используется для расчета ICS.

type icEvaluator.m

classdef icEvaluator
% icEvaluator Evaluates Initial Conditions data at requested locations
% ICE = icEvaluator(R) Creates an initial conditions evaluator from a 
% results object. The evaluator provides a function that can be called at 
% specific locations within the geometry. This class customized to represent
% a stationary solution.
%
% icEvaluator methods:
%    computeIC - Computes ICs at locations within the geometric domain
%

% Copyright 2015 The MathWorks, Inc.              
    properties
        Results;           
    end

    methods
        function obj = icEvaluator(R)                 
            obj.Results = R;                                                                                                                                 
        end
                
        function ic = computeIC(obj,locations)
            is2d = size(obj.Results.Mesh.Nodes, 1) == 2;
            if is2d
                querypts = [locations.x; locations.y];
            else
                  querypts = [locations.x; locations.y; locations.z];
            end
            numeqns = size(obj.Results.NodalSolution,2);
            if numeqns == 1 
                ic = interpolateSolution(obj.Results, querypts);
            else       
                ic = zeros(numeqns, numel(locations.x));
                for i = 1:numeqns
                    ic(i,:) = interpolateSolution(obj.Results, querypts, i);
                end             
            end           
        end       
    end    
end

Документация