MATLABで学ぶ振動工学 メタマテリアル・周期構造 3.0

 

注意:make_AR_AL.mには間違いがありますので、下記の縮約なしの記事を参考にしてください。

MATLABで学ぶ振動工学 メタマテリアル・周期構造 4.0(内部自由度の縮約なし)
D&D2019で発表された「周期構造の弾性波バンドギャップに基づくシェル構造の振動低減」を参考にメタマテリアルについて解説したいと思います。FEMモデルやMK行列など必要なものがすべて載せています。

 

はじめに

前回、1次元モデルで分散関係について説明しました。そして、「分散関係がメタマテリアルや周期構造物のバンドギャップ生成の根本的な考え方です」という説明をさせていただきました。

本記事では2次元モデルの説明をします。

いざ2次元モデルを説明しようとストーリーを考えて、がちがちの理論展開した記事を書こうとしたのですがやめました!おそらく、Blochの定理まででつまづいてしまうかなという気がしました。

ただ、Blochの定理までについては後日解説したいと思います。(いつになるかはわかりませんが…)

”メタマテリアルを勉強したい!でも、できてない!”という方は、勉強する取っ掛かりが欲しいのかなと思っています。私の役割はハードルをぐーーーっと下げて、自分の手でコードを実行してもらうことだと思いますので、ここでは難しい部分の説明は避けて、とりあえずコードを動かして結果が出力できるものを提供したいと思います。

ちなみに、この記事は下記文献を参考にしています。

2019年 機械学会(Dynamics and Design Conference)
“周期構造の弾性波バンドギャップに基づくシェル構造の振動低減”
豊田中央研究所 富田直, 西垣英一, 表竜二
https://www.jstage.jst.go.jp/article/jsmedmc/2019/0/2019_109/_article/-char/ja/

天下の豊田中研様の論文です。

しかし、この参考文献だけではプログラムを作ることはできません。論文にはよくあることですが、”論文として発表したいけど、内容は教えたくない。でも、すごいことしているように思われたい”ということが根底にあるかと思います。本参考文献もそんな感じです。

プログラムを書けるぐらい論文を紐解くには、関連する論文をざっと10冊程度読まないと厳しい、といのが実状かと思いますが、私は趣味で勉強しているので、ここでは余すことなく情報展開していきます。

メタマテリアルの構造・モデル

参考文献からこの記事用にモデルを作成しました。
なお、参考文献よりも節点(ノード)を大幅に削減して、自由度を低減しています。(下記リンクでモデルを提供するには容量を低減する必要があったため)

また、MK行列は自作のシェルプログラムで構築しているので、論文の結果とは一致しませんので、ご了承ください。

 

この記事で使用するユニットセルのモデルは下記です。このユニットセルが連続して並んでいる平板がメタマテリアルだと思ってください。ただし、メタマテリアルの性能はユニットセルで決まるので、下記モデルだけを扱います。

FEMモデルが欲しい方・確認したい方は下記を参照ください。

MATLABで学ぶ振動工学 メタマテリアル・周期構造 3.1 (ユニットセルのFEMモデル)
メタマテリアルの解説用に使用するFEMモデルです。

・重要
ちなみに、このユニットセルのノードID(節点番号)は下図のように外周部が1~24となりように割り振っています。(ここら辺の理由もちゃんと説明した方が理解してもらえそうですが、深夜になってしまったので本日は割愛します。解説してほしい方はコメント頂けると幸いです。)

バンドギャップの計算方法

メタマテリアルや周期構造は下記の定理・条件を前提としています。

  • 周期境界条件(ボルン=フォン・カルマン境界条件)
  • シュレディンガー方程式
  • シュレディンガー方程式の並進対称性
  • ハミルトニアン演算子
  • ブロッホの定理
    などなど

上記については後日詳細をまとめたいと思います。(たぶん)

 

ざっくりブロッホの定理までを説明すると、メタマテリアルや周期構造は”1つのユニットセルだけを考えればOKです”ってことです。

分散特性の導出方法

1つだけのユニットセル(単位ユニットセル)の運動方程式は下記で表せます。

$$
[K-ω^2 M]q =f \tag{1}
$$

Mは質量行列、Kは剛性行列、qは一般化変位、fは一般化内力です。
ユニットセルにブロッホの定理を適用するために、下記のように自由度を分類します。

\(q= \left[
\begin{array}{ccccccccc}
q_{LB}^T & q_{RB}^T & q_{LT}^T & q_{RT}^T & q_{L}^T & q_{B}^T & q_{R}^T & q_{T}^T & q_{I}^T
\end{array}
\right]^T\)

\(f= \left[
\begin{array}{ccccccccc}
f_{LB}^T & f_{RB}^T & f_{LT}^T & f_{RT}^T & f_{L}^T & f_{B}^T & f_{R}^T & f_{T}^T & f_{I}^T
\end{array}
\right]^T\)

さらに、上図のqI、fI以外を境界自由度、qI、fIを内部自由度に再分類して、下記の表現を用います。

\(q= \left[
\begin{array}{cc}
q_{Bd}^T & q_{I}^T
\end{array}
\right]^T\)

\(f= \left[
\begin{array}{cc}
f_{Bd}^T & 0
\end{array}
\right]^T\)

そうすると、境界自由度と内部自由度を分けて運動方程式を再構築することができます。

\[
(\left[ \begin{array}{cc} K_{Bd,Bd} & K_{Bd,I} \\ K_{I,Bd} & K_{I,I} \end{array} \right]
-ω^2 \left[ \begin{array}{cc} M_{Bd,Bd} & M_{Bd,I} \\ M_{I,Bd} & M_{I,I} \end{array} \right])
\left[ \begin{array}{c} q_{Bd} \\ q_{I} \end{array} \right]
=\left[ \begin{array}{c} f_{Bd} \\ 0 \end{array} \right]
\]

さらに、内部自由度をモード縮約して分散曲線を得る方法を用いて、下記のように表すことができます。(参考文献:Zhou, C., Laine, J., Ichchou, M., Zine, A., “Multi-scale modelling for two dimensional periodic structures using a combined mode/wave based approach”, Computers & Structures , Vol. 154 (2015) pp.145 162.)

\[
\left[ \begin{array}{c} q_{Bd} \\ q_{I} \end{array} \right]=
\left[ \begin{array}{cc} I & 0 \\ Ψ_{Bd} & Ψ_{C} \end{array} \right]
\left[ \begin{array}{c} q_{Bd} \\ P_{C} \end{array} \right]=
B \left[ \begin{array}{c} q_{Bd} \\ P_{C} \end{array} \right]
\]

ここで、Ψcはユニットセルの境界を完全拘束したときのモードベクトルであり、ΨBdは下記です。

$$
Ψ_{Bd}=-K_{II}^{-1}K_{I,Bd}
$$

$$
M_{reduc}=B^{T}MB, K_{reduc}=B^{T}KB
$$

参考までに載せますが、境界を完全拘束したモデルは下記です。

MATLABで学ぶ振動工学 メタマテリアル・周期構造 3.4 (ユニットセルのFEMモデル 周辺完全拘束)
メタマテリアルの解説用に使用するFEMモデルです。

そして、ブロッホの定理を適用すると一般化変位は下式で表すことができます。

\(
q_{Bd}=AR
\left[ \begin{array}{c} q_{LB}  \\ q_{L}  \\ q_{B} \end{array} \right]
=\left[ \begin{array}{ccc}
I & 0 & 0 \\
λ_xI & 0 & 0 \\
λ_yI & 0 & 0 \\
λ_xλ_yI & 0 & 0 \\
0 & I & 0  \\
0 & 0 & I  \\
0 & λ_xI & 0  \\
0 & 0 & λ_yI  \\
\end{array} \right]
\left[ \begin{array}{c} q_{LB}  \\ q_{L}  \\ q_{B} \end{array} \right]
\)

一般化内力は下記のように表すことができます。

\(
AL f_{Bd}=
\left[ \begin{array}{c} f_{LB}  \\ f_{L}  \\ f_{B} \end{array} \right]
\left[ \begin{array}{ccc}
I & 0 & 0  \\
λ_xI^-1 & 0 & 0  \\
λ_yI^-1 & 0 & 0  \\
λ_x^-1λ_y^-1I & 0 & 0  \\
0 & I & 0  \\
0 & 0 & I  \\
0 & λ_x^-1I & 0  \\
0 & 0 & λ_y^-1I  \\
\end{array} \right]^T
f_{Bd}=0
\)

ここで、λx、λyは波数kxおよび波数ky、ユニットセルの長さLxおよびLyを用いて下記となります。

$$
λ_x =e^{jk_x L_x}, λ_y =e^{jk_y L_y},
$$

さらに、質量行列と剛性行列の自由度を縮約することができます。

$$
MM=AL M_{reduc} AR, KK=AL M_{reduc} AR
$$

このMM行列とKK行列の固有値(共振周波数)が分散曲線となります。

MATLABでの検証結果

まず、検証結果を下記に示します。なんとなく参考文献に載っているような図になりました。下図の赤点線で示すように560~650Hzあたりにバンドギャップがありますね~。

参考文献ではこのユニットのパラメータースタディをしており、どんな感じにするとバンドギャップが拡大できるのかを考察しています。

もう少し詳しく書こうと思いましたが、かなり時間がかかってしまったので、質問のある方はコメントかツイッターでお願いします。
回答は別途記事を作成する形にしたいと思っています。(図なしで説明するの難しいので)

参考文献のパラメータースタディについても後日アップするかもしれません。

MATLABコード

実行ファイル

clear all
clc

[M]=get_m;
[K]=get_k;
% % LとRの節点数
nL=5; nR=5;
% % TとBの節点数
nT=5; nB=5;

nAll=nL+nR+nT+nB+4;



Kbd=K(1:6*nAll,1:6*nAll);
KIbd=K(6*nAll+1:end,1:6*nAll);
KII=K(6*nAll+1:end,6*nAll+1:end);

invKII=pinv(KII);
invKII(isnan(invKII))=0;
Psi_bd=invKII*KIbd;

Psi_bd=-invKII*KIbd;

B=[eye(size(Psi_bd,2))
    Psi_bd];

Kc=K; Mc=M;
Kc(1:6*nAll,:)=[]; Mc(1:6*nAll,:)=[];
Kc(:,1:6*nAll)=[]; Mc(:,1:6*nAll)=[];

[Psi_C,D]=eig(full(Kc),full(Mc));
Psi_C=Psi_C(:,1:50);

B=[B, [zeros(size(B,1)-size(Psi_C,1),size(Psi_C,2)) ; Psi_C]];

K_reduc=B.'*K*B;
M_reduc=B.'*M*B;

Lx=150; Ly=150;
ux=linspace(0,pi,10); uy=linspace(0,pi,10); 

fn=cell(length(ux),length(uy));
j=sqrt(-1);
for ii1=1:length(ux)
    for ii2=1:length(uy)
        lamdax=exp(-j*ux(ii1));
        lamday=exp(-j*uy(ii2));
        
        KK_re=K_reduc(1:nAll*6, 1:nAll*6);
        MM_re=M_reduc(1:nAll*6, 1:nAll*6);
                
        [AR,AL]=make_AR_AL(nL,nB,size(KK_re,1),lamdax,lamday);
        
        KK=AL*KK_re*AR;
        MM=AL*MM_re*AR;
        
        [V, D]=eig(KK,MM);
        fn{ii1,ii2}=sqrt(diag(D))/(2*pi);
    end
end
%%
figure
plot(0*ones(1,length(fn{1,1})),abs(fn{1,1}),'o')
text(0,0,'Node1')
hold on
for ii1=2:10
    plot((ii1-1)*ones(1,length(fn{1,1})),abs(fn{ii1,1}),'o')
end
text(9,0,'Node13')
for ii1=1:10
    plot((ii1+10-1)*ones(1,length(fn{1,1})),abs(fn{10,ii1}),'o')
end
text(19,0,'Node14')

for ii1=1:10
    plot((ii1-1+20)*ones(1,length(fn{1,1})),abs(fn{10+1-ii1,10+1-ii1}),'o')
end
text(29,0,'Node17')
plot([9 9],[0 1500],'k--')
plot([19 19],[0 1500],'k--')
plot([29 29],[0 1500],'k--')
plot([0 29],[563 563],'r--')
plot([0 29],[651 651],'r--')
ylim([0 1500])
xlabel('境界の位置')
ylabel('Freq. [Hz]')

 

functionファイル

function [AR,AL]=make_AR_AL(nL,nB,LengMatrix,lamdax,lamday)

AR=zeros(LengMatrix,6 + nL*6 + nB*6);
AR(1:6,1:6)=eye(6);
AR(6+1:6*2,1:6)=lamdax*eye(6);
AR(6*2+1:6*3,1:6)=lamday*eye(6);
AR(6*3+1:6*4,1:6)=lamdax*lamday*eye(6);

AR(:, 6+1:6+nL*6)=[AR(1:6*4,6+1:6+nL*6)
eye(nL*6)
zeros(nL*6)
lamdax*eye(nL*6)
zeros(nL*6)
];

AR(:,6+nL*6+1:end)=[AR(1:6*4,6+nL*6+1:end)
zeros(nB*6)
eye(nB*6)
zeros(nB*6)
lamday*eye(nB*6)
];

AL=zeros(LengMatrix,6 + nL*6 + nB*6);
AL(1:6,1:6)=eye(6);
AL(6+1:6*2,1:6)=lamdax^-1*eye(6);
AL(6*2+1:6*3,1:6)=lamday^-1*eye(6);
AL(6*3+1:6*4,1:6)=lamdax^-1*lamday^-1*eye(6);

AL(:, 6+1:6+nL*6)=[AL(1:6*4,6+1:6+nL*6)
eye(nL*6)
zeros(nL*6)
lamdax^-1*eye(nL*6)
zeros(nL*6)
];

AL(:,6+nL*6+1:end)=[AL(1:6*4,6+nL*6+1:end)
zeros(nB*6)
eye(nB*6)
zeros(nB*6)
lamday^-1*eye(nB*6)
];
AL=AL.';

end

M行列

MATLABで学ぶ振動工学 メタマテリアル・周期構造 3.2 (質量行列)
メタマテリアルの解説用のユニットセルのM行列です。

K行列

MATLABで学ぶ振動工学 メタマテリアル・周期構造 3.3 (剛性行列)
メタマテリアルの解説用のユニットセルのK行列です。

コメント

  1. konan より:

    コメント失礼します.
    詳しい解説を公開していただき,大変勉強になりました.

    境界条件についてわからない点があったので,2点質問させていただけないでしょうか.

    1.記事中の例題ではユニットセルの境界を完全拘束していますが,これは内部自由度を縮約するために必要な境界条件という理解で正しいでしょうか?それともブロッホの定理を適用する上で必須になる境界条件なのでしょうか?

    2.数百節点のユニットセルにおいても,縮約有り無しで計算時間が変わってくるのでしょうか?

    縮約操作とブロッホの定理の切り分けができておらず,的はずれな質問でしたら申し訳ありません.
    モードベクトルでqを表す所で詰まってしまったため,縮約無しで分散関係を計算してみたいと考えています.

    お手すきの際にご教授いただけると幸いです.

    • konanさま

      コメントありがとうございます。

      > 1の質問について
      完全拘束は内部自由度を縮約するために必要な境界条件です。
      (Zhou著の参考文献を適用するために必要)

      > 2の質問について
      数百節点であれば計算時間が変わらないと思います。
      (考え方としては、内部自由度をモード次数に縮退するので、内部自由度>>モード次数でないと計算時間は変わらないと思います。)

      > 縮約操作無しでの分散関係の計算方法
      私自身は論文をトレースすることに注力していたので、縮約なしの求め方を理解できていませんが、上から6番目の式のMK行列で固有値で求められる気がします。
      時間を見つけて、上から6番目の式のMK行列で固有値と縮約操作した固有値が一致しているかを検証してみたいと思います。

      また何かございましたら気兼ねなくコメント頂けると幸いです。

      • konan より:

        ご返信ありがとうございます.よく理解できました.
        本記事で紹介して頂いた2つの文献を購入したので,自分でも手を動かして深堀りしていきたいと思います.

        畑違いの分野から音振担当となり困り果ててた時に偶然このHPを拝見し,明快でわかりやすい記事でいつも勉強させていただいています.

        お言葉に甘えて,またご質問させていただくことがあるかもしれませんが,何卒よろしくお願いいたします.

        • konanさま

          縮約なしで計算できました。
          これから文書化して、公開しますので少々おまちください。

          縮約なしでの検証で、本記事のプログラムの間違いを発見しました。。。。
          本記事のmake_AR_AL.mには間違いがありますので、ご注意ください。