注意：make_AR_AL.mには間違いがありますので、下記の縮約なしの記事を参考にしてください。
はじめに
前回、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モデルが欲しい方・確認したい方は下記を参照ください。
・重要
ちなみに、このユニットセルのノード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
$$
参考までに載せますが、境界を完全拘束したモデルは下記です。
そして、ブロッホの定理を適用すると一般化変位は下式で表すことができます。
\(
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コード
好評につき、編集したプログラムを下記で販売してみました。
コメント
コメント失礼します．
詳しい解説を公開していただき，大変勉強になりました．
境界条件についてわからない点があったので，2点質問させていただけないでしょうか．
1.記事中の例題ではユニットセルの境界を完全拘束していますが，これは内部自由度を縮約するために必要な境界条件という理解で正しいでしょうか？それともブロッホの定理を適用する上で必須になる境界条件なのでしょうか？
2.数百節点のユニットセルにおいても，縮約有り無しで計算時間が変わってくるのでしょうか？
縮約操作とブロッホの定理の切り分けができておらず，的はずれな質問でしたら申し訳ありません．
モードベクトルでqを表す所で詰まってしまったため，縮約無しで分散関係を計算してみたいと考えています．
お手すきの際にご教授いただけると幸いです．
konanさま
コメントありがとうございます。
> 1の質問について
完全拘束は内部自由度を縮約するために必要な境界条件です。
（Zhou著の参考文献を適用するために必要）
> 2の質問について
数百節点であれば計算時間が変わらないと思います。
（考え方としては、内部自由度をモード次数に縮退するので、内部自由度＞＞モード次数でないと計算時間は変わらないと思います。）
> 縮約操作無しでの分散関係の計算方法
私自身は論文をトレースすることに注力していたので、縮約なしの求め方を理解できていませんが、上から6番目の式のMK行列で固有値で求められる気がします。
時間を見つけて、上から6番目の式のMK行列で固有値と縮約操作した固有値が一致しているかを検証してみたいと思います。
また何かございましたら気兼ねなくコメント頂けると幸いです。
ご返信ありがとうございます．よく理解できました．
本記事で紹介して頂いた2つの文献を購入したので，自分でも手を動かして深堀りしていきたいと思います．
畑違いの分野から音振担当となり困り果ててた時に偶然このHPを拝見し，明快でわかりやすい記事でいつも勉強させていただいています．
お言葉に甘えて，またご質問させていただくことがあるかもしれませんが，何卒よろしくお願いいたします．
konanさま
縮約なしで計算できました。
これから文書化して、公開しますので少々おまちください。
縮約なしでの検証で、本記事のプログラムの間違いを発見しました。。。。
本記事のmake_AR_AL.mには間違いがありますので、ご注意ください。
MATLABパイセン様
コメント失礼いたします。
分かりやすい解説をしていただき参考にさせて頂いております。
ユニットセルのノードID（節点番号）は下図のように外周部が1～24となりように割り振っていますが、なぜそのような振り方になるのかご教示お願いしたいです。
また、FEMソフトによってはこのような接点番号の振り方にならないときがあるのですが、そのような時は、ソフトから出力される剛性行列等を変換することで対応可能でしょうか？
以上、お手すきの際にご教授お願いできますでしょうか？
モデルを作成するソフトにもよりますが、一般的なソフトでは節点番号の昇順にM行列およびK行列が作成されます。
このM行列およびK行列からK_{Bd,Bd}とK_{I,Bd}とK_{I,I}に分解したいということになります。
この分解にする際に、Bdの節点が1～24とわかっていると、K_{Bd,Bd}とK_{I,Bd}とK_{I,I}が簡単に分解することができます。
> FEMソフトによってはこのような接点番号の振り方にならないときがあるのですが、そのような時は、ソフトから出力される剛性行列等を変換することで対応可能でしょうか？
ソフトから出力される剛性行列などを変換することで対応可能です。
MATLABパイセン様
先日は質問に答えて頂きありがとうございました。
ご教授頂いた通り、Matlabでプログラムを組んでみたところ追加で疑問点があったため、質問させていただきます。
お手すきの際、ご教授頂けないでしょうか。
以前、noteの方でモデルを購入させて頂き、Femapの方から剛性行列、質量行列を出力し、以前ご教授頂いた通りmatlabで節点を変換して分散曲線を出力させていただきました。
その際、縮約なしのモデルではほぼMATLABパイセン様の記事と同様の分散曲線が得られたのですが、縮約ありで試してみると、ところどころ途切れたような分散関係となってしまいました。このことについて、出力した剛性行列がMATLABパイセン様のものと異なるためだということは分かったのですが、どのようにすれば同様の行列を出力できるのかわからず詰まってしまいました。
そのため、独自のシェルプログラムから出力した際、行列をどのような処理をされているのか可能な限りご教授いただでないでしょうか。
私自身、学生で有限要素法については学び始めたばかりなので初歩的な質問でしたらすみません、、
コメントありがとうございます。
有限要素法やメタマテリアルなどを熱心に勉強されている学生さんに再度コメント頂けたということで大変嬉しいです。
私も有限要素法を10年以上前に勉強して、特にシェル要素が非常に難しかった記憶があります。シェル要素は厚板や薄板、ハイブリッドの理論など複数あるので、FEMの解析ソフトによっても剛性行列の出力結果が異なります。10年以上前に勉強して作成したコードなので、残念ながら詳細に説明できるほど覚えておらず、お力になれないかと思います。
noteで紹介しているプログラムは縮約なしだったかと思います。
noteに縮約ありのプログラムを載せていないのは、私自身が「縮約ありのコードが絶対に正しい」と言えるほど理解できていないためです。”縮約あり”の結果で試行錯誤して”縮約なし”と同様の結果になったので、自己研鑽の学習としては「この程度で良いか」と思いやめてしまったという経緯があります。
また、個人的な意見として、「そもそも”縮約あり”で解く必要性があるか？」という疑問があります。縮約して解くメリットは計算が早くなることだと理解していますが、ユニットセル程度のモデル規模であれば、計算自体に時間がかからないので、メリットがないように感じています。それなら、わざわざ複雑な縮退の式を理解する必要がないのでは、と考えている次第です。
コメント頂いた回答としては、「縮約ありで試してみると、ところどころ途切れたような分散関係となってしまいました。」という点の本質が、剛性行列に起因しないのではないかと考えています。
有限要素法の最も簡単な理論はソリッド要素です（四面体要素や六面体要素）。noteでダウンロードできるシェル要素のモデルと類似するようなFEMモデルをソリッド要素で作成してみてはいかがでしょうか？
そのソリッド要素で”縮約なし”の分散曲線を解くと、シェル要素の”縮約なし”の分散曲線とほとんど同じような分布になるかと思います。
節点のIDの対応に応じて、MATLABプログラムを修正する必要はあるかと思いますが、節点とK行列の要素位置の対応がご自身の中で深く理解できるかと思いますので、問題に本質に近づくのではないかと考えています。
kiu様の期待していた回答ではないかと思いますが、何かあればまたコメント頂けると幸いです。
MATLABパイセン様
早速の返信ありがとうございます。
貴重なご意見頂きとても参考になります。
特に、縮約についての論文もある程度読み込んだつもりになっていましたが、途切れたような分散関係の本質が剛性行列に起因しないという点参考にさせて頂き、コードをブラッシュアップできないか試行錯誤させて頂こうと思います。（出力した剛性行列がMATLABパイセン様のものと異なるという点について、出力した質量行列がMATLABパイセン様のものと形も数値もほぼ完全に一致していた点や、Nastranの設定で仮想回転剛性の設定を変えることでMATLABパイセン様の剛性行列に値が近づく点で、剛性行列が原因ではないかと疑っていました。）
私自身、元は主に形状の最適化の研究をしているため、計算時間の短縮を行える縮約ありの計算にメリットを感じています。そのため、縮約ありのプログラムを完成させたいので、おっしゃる通り、ソリッド要素でK行列についての理解を深め、様々なFEMを試しつつもう少し頑張ってみようと思います。