Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

Simscape に対するカンチレバー トラス構造の解析モデル

この例では、静的領域と周波数領域の両方で、通常のカンチレバー トラス構造におけるジョイントの変位に対するパラメーター化された解析式を求める方法を示します。静的なケースの解析式は厳密になります。周波数応答関数の表現は実際の周波数応答を近似して低次元化したものです。

この例では、次の Symbolic Math Toolbox™ 機能を使用します。

モデル パラメーターの定義

この例の目標は、カンチレバー トラス構造内の特定の桟で変位 d を一様な断面積パラメーター A に解析的に関連付けることです。支配方程式は Md2xdt2+Kx=F で表されます。ここで、d はカンチレバー右上のジョイントの負荷から発生します。トラスは左端のジョイントで壁に取り付けられています。

トラスの桟は、一様密度の線形弾性材料で作成されています。青で強調表示されている (図を参照) 桟の断面積 A は変化する可能性があります。他の桟の一様な断面積を含むその他すべてのパラメーターは、固定されています。微小な線形変位の仮定を使用して、先端の変位を取得します。

まず、トラスの固定パラメーターを定義します。

トラスの長さと高さおよび上部水平トラスの桟の数を指定します。

L = 1.5;
H = 0.25;
N = 3;

トラスの桟の材料の密度と弾性率を指定します。

rhoval = 1e2;
Eval = 1e7;

その他のトラスの桟の断面積を指定します。

AFixed = 1;

次に、特定のトラス要素の局所剛性行列を定義します。

局所剛性行列 k を計算し、シンボリック関数として表します。トラス要素の力と変位は局所剛性行列を介して関連付けられます。トラス要素の各ノードには、水平方向と垂直方向の変位の 2 つの自由度があります。各トラス要素の変位は [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) = 

(σ2σ1-σ2-σ1σ1σ3-σ1-σ3-σ2-σ1σ2σ1-σ1-σ3σ1σ3)where  σ1=AEsin(2t)2l  σ2=AEcos(t)2l  σ3=AEsin(t)2l

同様に、質量行列を計算します。

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 の列には 4 つの要素があり、以下のパラメーターを表します。

  • 長さ (l)

  • 横軸に対する角度 (t)

  • 開始ジョイントのインデックス

  • 終了ジョイントのインデックス

bars = zeros(2*N+2*(N-1),4);

上部と対角方向の桟を定義します。対象の桟は、(N+1) 番目の桟または番号 4 の桟で、左から 1 番目の対角方向の桟です。

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

桟を集めてシンボリック グローバル行列を作成

トラス構造には 7 つのノードがあります。各ノードには 2 つの自由度があります (水平方向と垂直方向の変位)。この系の自由度の合計数は 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 方程式の作成

A の関数として変位 d の解析解を求めてプロットします。ここで、K\F は全ジョイントの変位を表し、c は特定の変位を選択します。まず、簡潔にするために数値を 16 桁で表したシンボリック解を示します。

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

-0.00000001118033988749895504.7737197247586A2+384.7213595499958A+22.3606797749979A1.28A+0.8944271909999158

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

Figure contains an axes object. The axes object with xlabel Cross-Section, A, ylabel Displacement, d contains an object of type functionline.

結果を、Simscape ブロック内部で使用する Simscape 方程式 ssEq に変換します。得られた Simscape 方程式を表示するために、セミコロンを削除します。

syms d
ssEq = simscapeEquation(d,d_Static);

ssEq の式の最初の 120 文字を表示します。

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));

table でデータを可視化します。

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) を使用して計算を高速化します。

  • 伝達関数 H(s,A)=c(s2M(A)+K(A))-1F のパデ近似を、級数の最初の 3 つのモーメントを使用して、周波数 s = 0 の周辺で求めます。伝達関数が与えられると、パデ近似モーメントを c(-K(A)-1M(A))jK(A)-1F として計算できるということです。ここで、j{0,2,4,6,...} はべき級数項 {1,s2,s4,s6,...} に対応します。この例では、HApprox(s,A) を最初の 3 つの項の和として計算します。この方法は、伝達関数の次数を下げる非常に基本的な手法です。

  • さらに計算を高速化するには、AFixed の周辺で A のテイラー級数展開の各モーメント項の内項を近似します。

vpa を使用して計算を高速化します。

digits(32);
K = vpa(K);
M = vpa(M);

K の LU 分解を計算します。

[LK,UK] = lu(K);

補助変数と補助関数を作成します。

iK = @(x) UK\(LK\x);
iK_M = @(x) -iK(M*x);
momentPre = iK(F);

最初の 3 つのモーメントに対応する周波数級数を定義します。ここで、s は周波数です。

syms s
sPowers = [1 s^2 s^4];

前の計算の d_Static と同じ、1 番目のモーメントを設定します。

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;

最初の 3 つのモーメントを表示します。式が大きくなるため、一部分のみを表示することもできます。

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'                                                                                                                                                                                                                                                  }

周波数応答を As の数値を含む MATLAB 関数に変換します。Symbolic Math Toolbox なしで結果の 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');

Figure contains an axes object. The axes object with xlabel Frequency, ylabel Frequency Response contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Symbolic, Numeric.

選択した区間では、解析と数値の曲線が近接していますが、一般に、s = 0 の近傍以外の周波数では、曲線は離れたものになります。また、A の値が AFixed と大きく異なる場合、テイラー近似の誤差は大きくなります。