MATLABで学ぶ振動工学メタマテリアル・周期構造 3.0 - MATLABパイセンが教える振動・騒音・音響・機械工学

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

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

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

はじめに
メタマテリアルの構造・モデル
バンドギャップの計算方法
MATLABでの検証結果
MATLABコード
1. 関連

はじめに

前回、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モデルが欲しい方・確認したい方は下記を参照ください。

404 NOT FOUND | MATLABパイセンが教える振動・騒音・音響・機械工学

Mechanical Engineering University

・重要
ちなみに、このユニットセルのノード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
$$

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

404 NOT FOUND | MATLABパイセンが教える振動・騒音・音響・機械工学

Mechanical Engineering University

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

$
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コード

好評につき、編集したプログラムを下記で販売してみました。

MATLABで学ぶ振動工学　メタマテリアル・周期構造（平面構造・シェル構造）｜バネマスさん

１．はじめに最近は「メタマテリアル」というワードを良く聞きますよね。私自身もこの分野に非常に可能性を感じており、自己研鑽も兼ねて勉強しています。メタマテリアルは「周期構造」や「周期構造物」とも呼ばれ、分散関係を導出するとバンドギャップを持つ構造を指します。一般的に、振動工学や振動モード解析などの分野では、...