はじめに
前回の記事で初コメントを頂きました。ありがとうございます。
コメント
①縮約有り無しで計算時間が変わってくるのでしょうか？
②縮約無しで分散関係を計算してみたいです
①に関して、「計算時間は変わらないと思う」と回答してしまいましたが、おそらく変わります。文頭で記載したように、もともとのプログラムが間違っていたので、そこを修正すると、縮退した方が計算時間が早いというのが想像できます。ごめんなさい。
②に関しては、本記事で回答させていただきます。
メタマテリアルの構造・モデル
以前の記事の繰り返しになりますが、モデルについて紹介します。なお詳細な説明については、前回の記事（文頭のリンク）をご確認ください。
この記事で使用するユニットセルのモデルは下記です。このユニットセルが連続して並んでいる平板がメタマテリアルだと思ってください。ただし、メタマテリアルの性能はユニットセルで決まるので、下記モデルだけを扱います。
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]
\]
そして、ブロッホの定理を適用すると一般化変位は下式で表すことができます。
ここからが前回と違う
\[
(\left[ \begin{array}{c} q_{Bd} \\ q_{I} \end{array} \right]
=AR
\left[ \begin{array}{c} q_{LB} \\ q_{L} \\ q_{B} \\ q_{I} \end{array} \right]
=\left[ \begin{array}{cccc}
I & 0 & 0 & 0 \\
λ_xI & 0 & 0 & 0 \\
λ_yI & 0 & 0 & 0 \\
λ_xλ_yI & 0 & 0 \\
0 & I & 0 & 0 \\
0 & 0 & I & 0 \\
0 & λ_xI & 0 & 0 \\
0 & 0 & λ_yI & 0 \\
0 & 0 & 0 & I \\
\end{array} \right]
\left[ \begin{array}{c} q_{LB} \\ q_{L} \\ q_{B} \\ q_{I} \end{array} \right]
\]
一般化内力は下記のように表すことができます。
\(
AL f=
\left[ \begin{array}{cccc}
I & 0 & 0 & 0 \\
λ_xI^-1 & 0 & 0 & 0 \\
λ_yI^-1 & 0 & 0 & 0 \\
λ_x^-1λ_y^-1I & 0 & 0 & 0 \\
0 & I & 0 & 0 \\
0 & 0 & I & 0 \\
0 & λ_x^-1I & 0 & 0 \\
0 & 0 & λ_y^-1I & 0 \\
0 & 0 & 0 & I \\
\end{array} \right]^T
f=0
\)
ここで、λx、λyは波数kxおよび波数ky、ユニットセルの長さLxおよびLyを用いて下記となります。
$$
λ_x =e^{jk_x L_x}, λ_y =e^{jk_y L_y},
$$
さらに、質量行列と剛性行列の自由度を縮約することができます。
$$
MM=AL M AR, KK=AL M AR
$$
このMM行列とKK行列の固有値（共振周波数）が分散曲線となります。
MATLABでの検証結果
まず、検証結果を下記に示します。なんとなく参考文献に載っているような図になりました。下図の赤点線で示すように330～380Hzあたりと600～700Hzあたりにバンドギャップがありますね～。
参考文献ではこのユニットのパラメータースタディをしており、どんな感じにするとバンドギャップが拡大できるのかを考察しています。
MATLABコード
実行ファイル
好評につき、編集したプログラムを下記で販売してみました。
コメント