В этом примере показано, как найти параметризованные аналитические выражения для смещения соединения тривиальной конструкции фермы-консоли как в статической, так и в частотной областях. Аналитическое выражение для статического случая является точным. Выражение для функции частотной характеристики является приблизительной версией фактической частотной характеристики с уменьшенным порядком.
В этом примере используются следующие возможности символьного математического Toolbox™:
Создание кода Simscape с помощью simscapeEquation
Разложение логической единицы с использованием lu
Создание функции MATLAB с использованием matlabFunction
Цель этого примера состоит в том, чтобы аналитически связать смещение d с параметром A однородной площади поперечного сечения для конкретного стержня в конструкции фермы-консоли. Управляющее уравнение представлено = F. Здесь d является результатом нагрузки в верхнем правом соединении консоли. Ферма крепится к стене в крайних левых соединениях.
Стержни ферм выполнены из линейного эластичного материала с равномерной плотностью. Площадь поперечного сечения A полосы, выделенная синим цветом (см. рисунок), может изменяться. Все остальные параметры, включая однородные площади поперечного сечения других стержней, являются фиксированными. Для получения смещения наконечника используйте предположения о малом и линейном смещении.
Сначала определите фиксированные параметры фермы.
Укажите длину и высоту фермы, а также количество верхних горизонтальных стержней фермы.
L = 1.5; H = 0.25; N = 3;
Задайте плотность и модуль упругости материала стержня фермы.
rhoval = 1e2; Eval = 1e7;
Задайте площадь поперечного сечения других стержней фермы.
AFixed = 1;
Затем определите локальную матрицу жесткости конкретного элемента фермы.
Вычислите локальную матрицу жесткости k и выразите ее как символическую функцию. Силы и перемещения элемента фермы связаны через локальную матрицу жесткости. Каждый узел элемента фермы имеет две степени свободы - горизонтальное и вертикальное смещение. Смещение каждого элемента фермы представляет собой матрицу колонн, соответствующую [x_node1,y_node1,x_node2,y_node2]. Полученная матрица жесткости представляет собой матрицу 4 на 4.
syms A E l t real G = [cos(t) sin(t) 0 0; 0 0 cos(t) sin(t)]; kk = A*E/l*[1 -1;-1 1]; k = simplify(G'*kk*G); localStiffnessFn = symfun(k,[t,l,E])
localStiffnessFn(t, l, E) =
Вычислите массовую матрицу аналогичным образом.
syms rho
mm = A*rho*l/6*[2 1;1 2];
m = simplify(G'*mm*G);
localMassFn = symfun(m,[t,l,rho]);Теперь определите стержни фермы как матрицу bars. Индексы стержней, стартовых соединений и концевых соединений показаны на следующем рисунке.
Количество строк матрицы bars - количество стержней фермы. Столбцы bars имеют четыре элемента, представляющих следующие параметры:
Длина (l)
Угол относительно горизонтальной оси (t)
Индекс стартового соединения
Индекс конечного соединения
bars = zeros(2*N+2*(N-1),4);
Определите верхнюю и диагональную полосы. Заметим, что интересующей полосой является (N + 1) -я полоса или полоса номер 4, которая является первой диагональной полосой слева.
for n = 1:N lelem = L/N; bars(n,:) = [lelem,0,n,n+1]; lelem = sqrt((L/N)^2+H^2); bars(N+n,:) = [lelem,atan2(H,L/N),N+1+n,n+1]; end
Определите нижнюю и вертикальную планки.
for n = 1:N-1 lelem = L/N; bars(2*N+n,:) = [lelem,0,N+1+n,N+1+n+1]; lelem = H; bars(2*N+N-1+n,:) = [lelem,pi/2,N+1+n+1,n+1]; end
Ферменная конструкция имеет семь узлов. Каждый узел имеет две степени свободы (горизонтальное и вертикальное смещения). Общее число степеней свободы в этой системе равно 14.
numDofs = 2*2*(N+1)-2
numDofs = 14
Матрица массы системы M и матрица жесткости системы K - символьные матрицы. Соберите матрицы локальных элементов из каждой полосы в соответствующую глобальную матрицу.
K = sym(zeros(numDofs,numDofs)); M = sym(zeros(numDofs,numDofs)); for j = 1:size(bars,1) % Construct stiffness and mass matrices for bar j. lelem = bars(j,1); telem = bars(j,2); kelem = localStiffnessFn(telem,lelem,Eval); melem = localMassFn(telem,lelem,rhoval); n1 = bars(j,3); n2 = bars(j,4); % For bars that are not within the parameterized area A, set the stiffness % and mass matrices to numeric values. Note that although the values % kelem and melem do not have symbolic parameters, they are still % symbolic objects (symbolic numbers). if j ~= N+1 kelem = subs(kelem,A,AFixed); melem = subs(melem,A,AFixed); end % Arrange the indices. ix = [n1*2-1,n1*2,n2*2-1,n2*2]; % Embed local elements in system global matrices. K(ix,ix) = K(ix,ix) + kelem; M(ix,ix) = M(ix,ix) + melem; end
Устранение степеней свободы, присоединенных к стене.
wallDofs = [1,2,2*(N+1)+1,2*(N+1)+2]; % DoFs to eliminate
K(wallDofs,:) = [];
K(:,wallDofs) = [];
M(wallDofs,:) = [];
M(:,wallDofs) = [];F - вектор нагрузки с нагрузкой, указанной в отрицательном Y направление на самое верхнее правое соединение.
F = zeros(size(K,1),1); F(2*N) = -1;
Извлечение Y в крайнем верхнем правом соединении создайте вектор выбора.
c = zeros(1,size(K,1)); c(2*N) = 1;
Поиск и печать аналитического решения для смещения d как функция A. Здесь, K\F представляет смещения на всех соединениях и c выбирает конкретные смещения. Сначала отобразите символьное решение, представляющее числовые значения, используя 16 цифр для компактности.
d_Static = simplify(c*(K\F)); vpa(d_Static,16)
ans =
fplot(d_Static,[AFixed/10,10*AFixed]); hold on; xlabel('Cross-Section, A'); ylabel('Displacement, d'); title('')
Преобразование результата в уравнение Simscape ssEq для использования внутри блока Simscape. Чтобы отобразить результирующее уравнение Simscape, удалите точку с запятой.
syms d
ssEq = simscapeEquation(d,d_Static);Отображение первых 120 символов выражения для ssEq.
strcat(ssEq(1:120),'...')ans = 'd == (sqrt(5.0)*(A*2.2e+2+A*cos(9.272952180016122e-1)*2.0e+2+sqrt(5.0)*A^2*1.16e+2+sqrt(5.0)*1.0e+1+A^2*cos(9.2729521800...'
Сравнение точного аналитического решения и числового решения в диапазоне A значения параметров. Значения образуют последовательность из AFixed кому 5*AFixed с шагом 1.
AParamValues = AFixed*(1:5)'; d_NumericArray = zeros(size(AParamValues)); for k=1:numel(AParamValues) beginDim = (k-1)*size(K,1)+1; endDim = k*size(K,1); % Create a numeric stiffness matrix and extract the numeric solution. KNumeric = double(subs(K,A,AParamValues(k))); d_NumericArray(k) = c*(KNumeric\F); end
Вычисление символьных значений по AParamValues. Для этого вызовите subs функция, а затем преобразовать результат в числа с двойной точностью, используя double.
d_SymArray = double(subs(d_Static,A,AParamValues));
Визуализация данных в таблице.
T = table(AParamValues,d_SymArray,d_NumericArray)
T=5×3 table
AParamValues d_SymArray d_NumericArray
____________ ___________ ______________
1 -4.6885e-06 -4.6885e-06
2 -4.5488e-06 -4.5488e-06
3 -4.5022e-06 -4.5022e-06
4 -4.4789e-06 -4.4789e-06
5 -4.4649e-06 -4.4649e-06
Параметрическое представление для частотной характеристики H(s,A) может быть эффективным для быстрого изучения влияния параметра A не прибегая к потенциально дорогостоящим числовым вычислениям для каждого нового значения A. Однако найти точное решение в закрытой форме для большой системы с дополнительными параметрами часто невозможно. Поэтому обычно аппроксимируется решение для таких систем. Этот пример аппроксимирует параметрическую частотную характеристику вблизи нулевой частоты следующим образом:
Ускорение вычислений с помощью арифметики переменной точности (vpa).
Найти приближение Паде передаточной функции K (A)) -1F вокруг частотыs = 0 используя первые три момента серии. Идея заключается в том, что, учитывая передаточную функцию, можно вычислить аппроксимационные моменты Паде как (A) -1F, j∈{0,2,4,6,...} соответствуют степенных рядов HApprox(s,A) как сумма первых трех слагаемых. Этот подход является очень основным методом уменьшения порядка передаточной функции.
Чтобы еще больше ускорить вычисления, аппроксимируйте внутренний член каждого члена момента в терминах расширения серии Тейлора A вокруг AFixed.
Использовать vpa для ускорения вычислений.
digits(32); K = vpa(K); M = vpa(M);
Вычислить разложение LU K.
[LK,UK] = lu(K);
Создание вспомогательных переменных и функций.
iK = @(x) UK\(LK\x); iK_M = @(x) -iK(M*x); momentPre = iK(F);
Определите частотный ряд, соответствующий первым трем моментам. Здесь, s - частота.
syms s
sPowers = [1 s^2 s^4];Установите первый момент, который совпадает с d_Static в предыдущем вычислении.
moments = d_Static;
Вычислите оставшиеся моменты.
for p=2:numel(sPowers) momentPre = iK_M(momentPre); for pp=1:numel(momentPre) momentPre(pp) = taylor(momentPre(pp),A,AFixed,'Order',3); end moment = c*momentPre; moments = [moments;moment]; end
Объединение элементов момента для создания приблизительной аналитической частотной характеристики Happrox.
HApprox = sPowers*moments;
Отображение первых трех моментов. Поскольку выражения велики, их можно отобразить только частично.
momentString = arrayfun(@char,moments,'UniformOutput',false); freqString = arrayfun(@char,sPowers.','UniformOutput',false); table(freqString,momentString,'VariableNames',{'FreqPowers','Moments'})
ans=3×2 table
FreqPowers Moments
__________ _____________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________
{'1' } {'-(5^(1/2)*(220*A + 200*A*cos(8352332796509007/9007199254740992) + 116*5^(1/2)*A^2 + 10*5^(1/2) + 40*A^2*cos(8352332796509007/9007199254740992) + 40*A^2 + 20*5^(1/2)*A + 152*5^(1/2)*A^2*cos(8352332796509007/9007199254740992) + 36*5^(1/2)*A^2*cos(8352332796509007/4503599627370496)))/(200000000*A*(A - A*cos(8352332796509007/4503599627370496) - 5^(1/2)*cos(8352332796509007/9007199254740992) + 5^(1/2)))'}
{'s^2'} {'0.000000000084654421099019119162064316556812*(A - 1)^2 - 0.000000000079001242991597407593795324876303*A + 0.0000000004452470882909076005654298481594' }
{'s^4'} {'0.000000000000012982169746567833536274593916742*A - 0.000000000000015007141035503798552656353179406*(A - 1)^2 - 0.000000000000045855285883001825767717087472451' }
Преобразование частотной характеристики в функцию MATLAB, содержащую числовые значения для A и s. Результирующую функцию MATLAB можно использовать без панели символьных математических инструментов.
HApproxFun = matlabFunction(HApprox,'vars',[A,s]);Изменение частоты от 0 кому 1 в logspace, а затем преобразовать диапазон в радианы.
freq = 2*pi*logspace(0,1);
Вычислить числовые значения частотной характеристики для A = AFixed*perturbFactor, то есть для небольшого возмущения вокруг A.
perturbFactor = 1 + 0.25; KFixed = double(subs(K,A,AFixed*perturbFactor)); MFixed = double(subs(M,A,AFixed*perturbFactor)); H_Numeric = zeros(size(freq)); for k=1:numel(freq) sVal = 1i*freq(k); H_Numeric(k) = c*((sVal^2*MFixed + KFixed)\F); end
Вычислите символьные значения частотной характеристики в диапазоне частот и A = perturbFactor*AFixed.
H_Symbolic = HApproxFun(AFixed*perturbFactor,freq*1i);
Постройте график результатов.
figure loglog(freq/(2*pi),abs(H_Symbolic),freq/(2*pi),abs(H_Numeric),'k*'); xlabel('Frequency'); ylabel('Frequency Response'); legend('Symbolic','Numeric');
Аналитические и числовые кривые близки в выбранном интервале, но в целом кривые расходятся для частот за пределами окрестности s = 0. Кроме того, для значений A далеко от AFixed, аппроксимация Тейлора вносит большие ошибки.