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

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

Цель этого примера состоит в том, чтобы аналитически связать смещение, d, к универсальному параметру площади поперечного сечения, A, для конкретной панели в консольной структуре связки. Здесь, d следует из загрузки в правильном верхнем соединении консоли. Связка присоединена к стене слева большинство соединений.

Панели связки характеризуются линейным эластичным материалом с универсальной плотностью. Кроме A, фиксируются все другие параметры включая площади поперечного сечения (универсальная форма) других панелей. Смещение совета получено с помощью небольших (линейных) предположений смещения.

Во-первых, задайте следующие фиксированные параметры связки.

Длина и высота связки и количество главной горизонтали связывают панели:

L = 1;
H = 0.25;
N = 2;

Плотность и модуль эластичности материала панели связки:

rhoval = 1e2;
Eval = 1e7;

Область поперечного сечения панели связки:

AFixed = 1;

Затем задайте локальную матрицу жесткости плоского элемента связки.

Вычислите локальную матрицу жесткости, k, и выразите ее как символьную функцию. Столбцы матрицы жесткости соответствуют [x_node1, y_node1, x_node2, y_node2];

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{где}\\ \\ \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}{\mathrm{потому\; что}\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 = zeros(2*N+2*(N-1),4);

Задайте верхние и двойные косые черты. Обратите внимание на то, что панель интереса является (N+1)th панелью, которая является первой двойной косой чертой слева.

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

Соберите панели в символьные глобальные матрицы

Количество степеней свободы (степень свободы) в этой системе

numDofs = 2*2*(N+1)-2

numDofs = 10

Системная большая матрица, 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 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 indices.
    ix = [n1*2-1,2*n1,n2*2-1,n2*2];
    % Embed local elements into system 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 из точного символьного решения для статического случая

Постройте аналитическое решение. Здесь, K\F представляет смещения во всех соединениях, и c выбирает конкретные смещения. Во-первых, распечатайте символьное решение, представляющее числовые значения с помощью 16 цифр для компактности.

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

ans = 
$$-\frac{0.000000025\hspace{0.17em}\left(77.03083505599866\hspace{0.17em}{A}^2+121.9646997739904\hspace{0.17em}A+20.0\right)}{A\hspace{0.17em}\left(1.28\hspace{0.17em}A+1.788854381999832\right)}$$

fplot(d_Static,[AFixed/10,10*AFixed]);
hold on;
xlabel('Cross-section, A');
ylabel('Displacement');
title('')

Преобразуйте результат в уравнение Simscape, ssEq, чтобы использовать в блоке Simscape. Чтобы отобразить получившееся уравнение Simscape, удалите точку с запятой.

syms d
ssEq = simscapeEquation(d,d_Static)

ssEq = 
'd == ((A*2.0e+1+sqrt(5.0)*A^2*4.0+A^2*cos(9.272952180016122e-1)*4.8e+1+A^2*cos(1.854590436003224)*1.1e+1+A^2*3.7e+1+sqrt(5.0)*A*3.0e+1+sqrt(5.0)*A*cos(9.272952180016122e-1)*2.6e+1+sqrt(5.0)*A^2*cos(9.272952180016122e-1)*4.0+2.0e+1)*(-2.5e-8))/(A*(A-A*cos(1.854590436003224)-sqrt(5.0)*cos(9.272952180016122e-1)*2.0+sqrt(5.0)*2.0));'

Сравните числовые и символьные статические решения

Сравните точное аналитическое решение и числовое решение в области значений значений параметров 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

Вычислите символьные значения по Aparamrange. Для этого вызовите функцию 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           -1.784e-06      -1.784e-06  
         2          -1.6443e-06     -1.6443e-06  
         3          -1.5977e-06     -1.5977e-06  
         4          -1.5744e-06     -1.5744e-06  
         5          -1.5604e-06     -1.5604e-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'         '-(20*A + 4*5^(1/2)*A^2 + 48*A^2*cos(8352332796509007/9007199254740992) + 11*A^2*cos(8352332796509007/4503599627370496) + 37*A^2 + 30*5^(1/2)*A + 26*5^(1/2)*A*cos(8352332796509007/9007199254740992) + 4*5^(1/2)*A^2*cos(8352332796509007/9007199254740992) + 20)/(40000000*A*(A - A*cos(8352332796509007/4503599627370496) - 2*5^(1/2)*cos(8352332796509007/9007199254740992) + 2*5^(1/2)))'
      's^2'       '0.00000000001931549865390535481955339404344*(A - 1)^2 - 0.000000000016488909600194499035418898203187*A + 0.000000000045018038433136039672236177596169'                                                                                                                                                                                                                                       
      's^4'       '0.0000000000000005984738130764688358429817612765*A - 0.00000000000000082093440076592193426384675490731*(A - 1)^2 - 0.0000000000000011768131579530973008654036198437'                                                                                                                                                                                                                         

Преобразуйте частотную характеристику на функцию 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, приближение Тейлора вводит большие ошибки.