Задайте параметры модели

Цель этого примера состоит в том, чтобы аналитически связать перемещение d с параметром площади равномерного поперечного сечения A для конкретного стержня в конструкции консольной фермы. Управляющее уравнение представлено $$M\frac{d^{2}x}{d{t}^2}+Kx=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) = 
$$\begin{array}{l}\left(\begin{array}{cccc}{\sigma }_2& {\sigma }_1& -{\sigma }_2& -{\sigma }_1\\ {\sigma }_1& {\sigma }_3& -{\sigma }_1& -{\sigma }_3\\ -{\sigma }_2& -{\sigma }_1& {\sigma }_2& {\sigma }_1\\ -{\sigma }_1& -{\sigma }_3& {\sigma }_1& {\sigma }_3\end{array}\right)\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\frac{A\hspace{0.17em}\text{E}\hspace{0.17em}\sin \left(2\hspace{0.17em}t\right)}{2\hspace{0.17em}l}\\ \\ \mathrm{ }{\sigma }_2=\frac{A\hspace{0.17em}\text{E}\hspace{0.17em}{\cos \left(t\right)}^2}{l}\\ \\ \mathrm{ }{\sigma }_3=\frac{A\hspace{0.17em}\text{E}\hspace{0.17em}{\sin \left(t\right)}^2}{l}\end{array}$$

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

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 имеют четыре элемента, которые представляют следующие параметры:

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;

Создайте Simscape уравнения из точного символьного решения для статического случая

Найдите и постройте график аналитического решения для d перемещения как функцию A. Здесь, K\F представляет смещения во всех соединениях и c выбирает конкретные перемещения. Сначала покажите символическое решение, представляющее числовые значения, используя 16 цифр для компактности.

d_Static = simplify(c*(K\F));
vpa(d_Static,16)

ans = 
$$-\frac{0.00000001118033988749895\hspace{0.17em}\left(504.7737197247586\hspace{0.17em}{A}^2+384.7213595499958\hspace{0.17em}A+22.3606797749979\right)}{A\hspace{0.17em}\left(1.28\hspace{0.17em}A+0.8944271909999158\right)}$$

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 для ускорения вычислений.

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 без Symbolic Math Toolbox.

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, приближение Тейлора вводит большие ошибки.