ドキュメンテーション

このページは前リリースの情報です。該当の英語のページはこのリリースで削除されています。

Simscape 解析モデルの作成

この例では、通常のカンチレバー トラス構造におけるジョイントの変位を静的領域と周波数領域でモデル化して Simscape で使用するために、パラメーター化された解析式を求めます。

静的なケースの解析式は厳密になります。周波数応答関数の表現は実際の周波数応答を近似して低次元化したものです。

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

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

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

トラスの桟は、一様密度の線形弾性材料の特徴を持ちます。他の桟の断面積 (一様) を含むその他すべてのパラメーターは、A を除き、固定されています。先端の変位は微小な (線形) 変位の仮定によって得られます。

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

トラスの長さと高さおよび上部の水平トラスの桟の数:

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

トラスの桟の材料の密度と弾性率:

rhoval = 1e2;
Eval = 1e7;

トラスの桟の断面積:

AFixed = 1;

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

局所剛性行列 (local stiffness matrix) 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) = 

(σ2σ1-σ2-σ1σ1σ3-σ1-σ3-σ2-σ1σ2σ1-σ1-σ3σ1σ3)ここで  σ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 の各列は次のとおりです。

  • 長さ (l)

  • 水平軸に対する角度 (t)

  • 開始ジョイント

  • 終了ジョイント

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

上部と対角方向の桟を定義します。対象の桟は、(N+1) 番目の桟で、左から 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

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

この系の自由度 (DoF) の数は以下です。

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 = 

-0.00000002577.03083505599866A2+121.9646997739904A+20.0A1.28A+1.788854381999832

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

結果を、Simscape ブロック内部で使用する Simscape 方程式 ssEq に変換します。得られた 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 の新しい値ごとに、演算負荷が高くなるような数値計算を利用しません。ただし、大きな系の厳密閉形式解を追加のパラメーターを使用して求めることは、多くの場合、不可能です。したがって、このような系については通常、解を近似します。この例では、次のようにゼロ周波数近くのパラメトリック周波数応答を近似します。

  • 伝達関数 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 つの項の和として計算します。この方法は、伝達関数の次数を下げる場合に適用される非常に基本的な手法です。文献で説明されているその他のよりロバストな手法は、この例の範囲外です。

  • 可変精度の演算 (vpa) を使用して計算を高速化します。

  • さらに計算を高速化するには、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];

1 番目のモーメントを設定します。これは、前に計算した 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'                                                                                                                                                                                                                         

周波数応答を As の数値を含む MATLAB 関数に変換します。結果の 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 と大きく異なる場合、テイラー近似の誤差は大きくなります。