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

Цель этого примера состоит в том, чтобы аналитически связать смещение d с универсальным параметром площади поперечного сечения для конкретной панели в консольной структуре связки. Управляющее уравнение представлено $$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)th панелью или панелью номер 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, приближение Тейлора вводит большие ошибки.