Transactions of JSCES, Paper No.20160024
熱変形制御を目的とする複合板のマルチスケールトポロジー最適化
Multiscale topology optimization for composite plates undergoing thermal displacement
西 紳之介
1，寺田 賢二郎
2，加藤 準治
3, 西脇 眞二4, 泉井 一浩4
Shinnosuke NISHI, Kenjiro TERADA, Junji KATO, Shinji NISHIWAKI and Kazuhiro IZUI
Shinnosuke NISHI, Kenjiro TERADA, Junji KATO, Shinji NISHIWAKI and Kazuhiro IZUI
1東北 大 学 大 学 院 工学 研究 科 土木工 学 専攻(〒980-0845仙台 市 青 葉区 荒 巻字 青 葉468-1)
2東北 大 学 災 害 科 学国 際研 究所(〒980-0845仙 台市 青葉 区 荒 巻字 青 葉468-1)
3東北大 学 大 学 院 工学研 究科 土 木工学 専 攻(〒980-8579仙台市 青 葉区 荒 巻 字 青葉6-6-06)
4京都大 学 大 学 院 工学研 究科 機 械理工 学 専攻(〒615-8540京 都 市 西京 区 京都 大 学 桂)
The present study proposes a multiscale topology optimization method of cross-section structures that control the macroscopic thermal deformations of a thick composite plates. The proposed optimization method is based on the two-scale composite plate model, in which thick plate theory is employed at macro- scale, while three-dimensional solids are assumed for in-plane periodic microstructures (unit cells) whose scale is comparable with the plate’s thickness. The process of the in-plane homogenization in a two-scale analysis corresponds to numerical plate testings (NPTs) to compute the macroscopic plate stiffness of the thick plate. In addition, the co-rotational formulation is employed to account for large displacement and ro- tation of the macroscopic plate structure, whereas material models used in in-plane unit cells, are assumed to be linearly elastic. The optimization problem that controls a nodal displacement of the macrostructure is se- lected as numerical examples, and their solutions are provided to demonstrate the capability of the proposed method.
Key Words: Multiscale topology optimization，Composite plates，Co-rotational formulation，In-plane homogenization, Thermal deformation
1. 緒 言
複合材料の性能は，用いる構成材料の材料特性はも ち ろ ん の こ と ，そ れ ら の 幾 何 形 状 や 組 み 合 わ せ 方 に 大 き く 依 存 す る た め ，求 め る 性 能 に 応 じ て 適 切 な ミ ク ロ 構 造 の 形 状 ，材 料 配 置 を 決 定 す る 必 要 が あ る ．求 め ら れ る 性 能 と し て は ，高 剛 性 化 を 図 る な ど し て 機 械 的 な 特 性 を 高 め る も の が 多 い が ，近 年 で は 複 数 の 機 能 を 持 た せ る 多 機 能 化 の 試 み が 広 がって い る[1–4]．そ の1例 と し て ア ク チュエ ー タ 機 能 が あ り，そ れ ら は 外 荷 重 の み な ら ず 温 度 や 電 流 な ど と いった 外 界 か ら の 作 用 に 対 し て 最 適 な 変 形 を 促 す と いった も の で あ る ．
そのような複合材料は材料の組み合わせ方が多様で あ る 一 方 ，そ の 加 工 技 術 の 限 界 も あ り，実 験 に よ る 設 計手 法では最適形状を発見しにく いという課題が指摘 さ れ て い る ．そ の た め 数 値 解 析 に よ る 構 造 設 計 に 関 す る 研 究 が 数 多 く な さ れ て お り，な か で も 複 雑 な 構 造 ，
∗ 原 稿 受 付2016年4月7日,改 訂 年 月 日2016年10月14日, 発 行 年 月 日2016年11月15日, c⃝2016年 日 本 計 算 工 学 会 ． Manuscript received, April 07, 2016; final revision, October 14, 2016;
published, November 15, 2016. Copyright c⃝2016 by the Japan Society for Computational Engineering and Science.
材 料 配 置 が 可 能 な 自 由 度 の 高 い 設 計 手 法 で あ る ト ポ ロジー最適化はこれまで多 くの研究例が報告され てき
た[5–8]．トポロジー最適化の研究はマクロ構造の設計
の み な ら ず，ミ ク ロ 構 造 を も 設 計 対 象 と す る 均 質 化 法 に基づくマルチスケールトポロジ ー最適化手法が報告 されるようになった[9–11]．しかし，これらの手法の多 くは，2次元あるいは3次元連続体モデルとして定式化 され，ミクロ構造・マクロ構造ともにソリッド要素を用 いた有限要素解析を前提としてい る場合がほとんどで あ り，面 内 方 向 に し か ミ ク ロ 構 造 が 周 期 的 に 配 置 さ れ な い 複 合 板 へ の 適 用 例 は 未 だ 報 告 さ れ て い な い ．こ の 理由として，Fig.1に示すように面内にのみ周期性を有 す る 複 合 板 の 場 合 ，前 述 の マ ル チ ス ケ ー ル 解 析 が 理 論 的拠り所としている通常の 均質化法の適用対象外 であ る こ と が 挙 げ ら れ る ．こ の よ う な 板 状 構 造 物 に 対 し て は，ミクロには3次元固体力学理論，マクロには厚板理 論を適用したマルチ スケール解析手法の 定式化が必要 とされる．そのような背景のもと，Teradaら[12, 13]は，
数学的均質化法のように板 厚や面内周期性の代表 寸法 を ゼ ロ に す る よ う な 操 作 を 行 わ ず，ユ ニット セ ル を 数 値供試体とみたてて数値試 験的に均質化剛性を算 出す る ，新 し い マ ル チ ス ケ ー ル 解 析 手 法 を 確 立 し た ．こ の
In-plane unit cell
: Representative length in in-plane domain
L : Representative length in in-plane domain
h
h x3=z3
x1 x2
y2 y3
y1
L >> l L
l l
Fig. 1 Composite plate with in-plane periodicity
手法により，板厚方向の材料の違いにより生じる，面外 せ ん 断 変 形 を 含 め た 複 雑 な 変 形 も 再 現 可 能 と な る ．ま た こ の 手 法 は 分 離 型 マ ル チ ス ケ ー ル 解 析 手 法[14]の 一 つ で あ る た め ，こ れ を ベ ー ス と す る 最 適 化 手 法[15, 16]
の ア ル ゴ リ ズ ム を 採 用 す る こ と が で き る ．
先述の複合板のマルチスケール解析[13]は微小変位・
微小回転の理論の枠組みで定式化されている．しかし，
ア ク チュエ ー タ な ど の 多 機 能 化 を 目 的 と し た 場 合 ，大 変 位・大 回 転 を 許 容 す る 幾 何 学 的 非 線 形 解 析 へ の 拡 張 は 必 須 で あ る ．こ の た め の 運 動 方 程 式 の 記 述 法 に は ， トータルLagrangian（TL），更新Lagrangian（UL），そし て共回転理論（co-rotational以後CR）がある．TL，ULを 用 い た 解 析 は ，そ れ ぞ れ 変 形 前 ，変 形 後 を 基 準 に つ り 合 い 方 程 式 を 解 く 手 法 で あ り，非 線 形 構 造 解 析 を 行 う 上 で 一 般 的 な 手 法 で あ る と い え る ．し か し ，こ れ ら の 手法 では有限変形理論に基づき非線形ひずみを用いて 定 式 化 が 行 わ れ て い る た め ，微 小 変 形 理 論 で 定 式 化 さ れて いる複合板のマルチスケール解析への適用が非常 に困難である．一方，CRによる定式化は，物体の運動 を剛 体運動と弾性変形による部分に区別して考える手 法 で あ り，材 料 構 成 則 に つ い て は 微 小 変 形 理 論 を 用 い る た め 複 合 板 の マ ル チ ス ケ ー ル 解 析 手 法 に よって 得 ら れ た 均 質 化 板 剛 性 を そ の ま ま 適 用 可 能 で あ る ．そ の た め ，TL，ULと 比 較 し て も 容 易 に 複 合 板 の 有 限 変 位 マ ル チ ス ケ ー ル 解 析 を 行 う こ と が 可 能 で あ る ．ま た ，共 回 転 理 論 に よ る 幾 何 学 的 非 線 形 解 析 は ，本 解 析 で 採 用 し た 板・シェル 構 造 解 析 へ の 適 用 性 に つ い て 多 く の 報
告例[17–19]があり，得られる解についても高い信頼性
を 有 し て い る こ と が 分 かって い る ．
そ こ で ，本 研 究 で は 面 内 に 周 期 的 な 非 均 質 性 を 有 す る 複 合 板 を 対 象 と し ，大 変 位・大 回 転 を 伴 う マ ク ロ 構 造 の 変 形 を 制 御 す る よ う な ミ ク ロ 構 造（ 面 内 ユ ニット セ ル ）の ト ポ ロ ジ ー を 決 定 す る マ ル チ ス ケ ー ル 最 適 化 計 算 手 法 の 開 発 を 目 的 と す る ．本 手 法 で は ，複 合 板 の マ ル チ ス ケ ー ル 解 析[13]を 採 用 し て ，面 内 ユ ニット セ ル に つ い て は3次 元 連 続 体 と し て モ デ ル 化 し て ソ リッ
N1
x1
M1
Neutral plane
1
u1
x3 1
x3
x3
u3(x
1,x
2,x
3)
=u
3(x
1,x
2) u3
x1
V1
x1 x3
In-plane unit cell
Fig. 2 Displacement field of thick plate theory
ド 要 素 を 用 い ，マ ク ロ 構 造 に つ い て は 厚 板 曲 げ 理 論 を 用 い た フ ラット シェル 要 素[20, 21]を 用 い る ．マ ク ロ 板 構 造 に つ い て は ，材 料 は 微 小 変 形 で あ る が ，大 変 位・
大回転問題を伴う問題にも 対応できるように共回 転理 論 を 採 用 す る ．最 適 化 計 算 に つ い て は ，分 離 型 の マ ル チ ス ケ ー ル 最 適 化 計 算 手 法[16]を 用 い ，温 度 変 化 が 与 えられたときにマクロ構造 の節点変位を指定した 方向 に 最 大 化 す る よ う な 最 適 化 問 題 を 設 定 す る ．ま た ，本 最適化計算では感度解析に 基づく最適化アルゴリ ズム を 採 用 す る ．こ の 際 ，感 度 に つ い て は 計 算 コ ス ト の 面 か ら 解 析 的 感 度 の 導 出 を 行 う．最 後 に 本 提 案 手 法 を 用 い て 得 ら れ た 数 値 解 析 例 を 提 示 し ，本 手 法 の 有 用 性 を 示 す．
2. 複合板のための面内均質化法
本 章 で は ，面 内 周 期 性 を 有 す る 複 合 板 を 対 象 と し た マ ル チ ス ケ ー ル 解 析 手 法 に つ い て 解 説 す る ．対 応 す る 2変 数 境 界 値 問 題 の 設 定 に 際 し て ，マ ク ロ 構 造 に は 厚 板 理 論 を ，ミ ク ロ 構 造（ 面 内 ユ ニット セ ル ）に つ い て は3次 元 固 体 力 学 を 適 用 し ，マ ク ロ 構 造 解 析 を 行 う た め の ，面 内 ユ ニット セ ル を 反 映 し た 均 質 化 板 剛 性 を 定 義 す る ．
ま ず，マ ク ロ 構 造 に 対 し て 厚 板 理 論 を 採 用 す る と ， マクロ座標 系x上に 働くマクロ 構造の変位場uはFig.2 を 参 照 し て 以 下 の よ う に な る ．
u1(x1,x2,x3)=¯u1(x1,x2)−x3ϕ1(x1,x2) u2(x1,x2,x3)=¯u2(x1,x2)−x3ϕ2(x1,x2) u3(x1,x2)=¯u3(x1,x2)
(1)
こ こ で ，¯uiは 中 心 面 上 の 変 位 で あ り，ϕiは 断 面 の 回 転 角を表している．また，図中のN，M，Vはそれぞれ板 の 断 面 に 働 く 垂 直 合 応 力 ，曲 げ モ ー メ ン ト 応 力 ，面 外 せん断合応力である．式(1)より，マクロひずみEは以 下 の よ う に 表 す こ と が で き る ．
y2
y3
y1
Mode1 Mode2 Mode3
Mode4 Mode5 Mode6
Mode7 Mode8
Fig. 3 Deformation modes imposed on unit cell with each macroscopic generalized strain
E11
E22
E33
2E12
2E23
2E31
=
∂¯u1
∂x1
−x3
∂ϕ1
∂x1
∂¯u2
∂x2
−x3
∂ϕ2
∂x2
0
∂¯u1
∂x2 +∂¯u2
∂x1 −x3
(∂ϕ2
∂x1 +∂ϕ1
∂x2
)
∂¯u3
∂x2
−ϕ2
∂¯u3
∂x1
−ϕ1
=
E˜1+x3E˜4 E˜2+x3E˜5
0 E˜3+x3E˜6
E˜7 E˜8
(2) こ こ で ，式 中 のE˜は マ ク ロ 一 般 化 ひ ず み で あ り，
E˜ ={
E˜1 E˜2 E˜3 E˜4 E˜5 E˜6 E˜7 E˜8}T
= {
∂¯u1
∂x1
∂¯u2
∂x2
∂¯u2
∂x1
+∂¯u1
∂x2
−∂ϕ1
∂x1
−∂ϕ2
∂x2
−
(∂ϕ2
∂x1
+∂ϕ1
∂x2
) ∂¯u3
∂x2
−ϕ2
∂¯u3
∂x1
−ϕ1
}T
(3) とおいた．これらの各マクロ一般化ひずみの変形図（マ ク ロ 変 形 モ ー ド ）をFig.3に 示 す．
次 に ，ミ ク ロ ひ ず みεは ，次 式 の よ う に 表 す．
ε=˜z ˜E+∂yu∗ (4) ここで，∂yu∗(y)は非均質に起因するミクロ擾乱ひずみ，
˜zは マ ク ロ 一 般 化 ひ ず み と ミ ク ロ ひ ず み を 結 び つ け る 行 列 で あ り，そ れ ぞ れ 次 式 の よ う に 定 義 し た ．
∂yu∗=
∂
∂y1
0 0
0 ∂
∂y2
0
0 0 ∂
∂y3
∂
∂y2
∂
∂y1
0
0 ∂
∂y3
∂
∂y2
∂
∂y3
0 ∂
∂y1
u∗1 u∗2 u∗3
(5)
˜z=
1 0 0 z3 0 0 0 0
0 1 0 0 z3 0 0 0
0 0 0 0 0 0 0 0
0 0 1 0 0 z3 0 0
0 0 0 0 0 −z1/2 1 0
0 0 0 0 0 −z2/2 0 1
(6)
ここで，u∗はy1,y2,y3の関数であり，面内ユニットセル 内に周期的に分布するミクロ擾乱変位である．˜zの(5,6) 成 分 ，(6,6)成 分 に ，−z1/2，−z2/2が 付 加 さ れ て い る が ， こ れ ら は ユ ニット セ ル の マ ク ロ ね じ り 変 形 を 再 現 す る た め に 導 入 さ れ た 成 分 で あ る[12, 13]．こ の と きyiとzi
と い う 二 つ の ミ ク ロ レ ベ ル の 座 標 系 を 導 入 し た が ，yi
は 均 質 化 の 過 程 で 用 い る も の で あ り，zi は 板 曲 げ 理 論 で 用 い る 運 動 学 的 な パ ラ メ ー タ で あ る ．ま た ，マ ク ロ ひ ず みEと ミ ク ロ ひ ず みεは 次 式 で 関 係 づ け ら れ て い る[12]．
E= 1 l1l2h
∫l1/2
−l1/2
∫ l2/2
−l2/2
∫h/2
−h/2εdy1dy2dy3
= 1 l1l2h
∫l1/2
−l1/2
∫ l2/2
−l2/2
∫h/2
−h/2
(˜z|(z1,z2)→(y1,y2)E˜+∂yu∗)
dy1dy2dy3
=z ˜E (7)
以上のマクロおよびミクロひずみの定義を用いると，
面 内 ユ ニット セ ル の 支 配 方 程 式 は 次 式 の よ う に な る ．
∂Tyσ=0 σ=C(
ε−εth) ε=˜z ˜E+∂yu∗=∂yw εth=α∆Tψ;ψ={
1 1 1 0 0 0}T
u∗: In-plane Y-periodic
(8)
こ こ で ，σ(y)は ミ ク ロ 応 力 ，wは ミ ク ロ 変 位 ，Cは 構 成材料の弾性係数行列であり，これらはすべてy1,y2,y3
の関数である．また，εthはミクロ熱膨張ひずみであり，
そ の 大 き さ は 熱 膨 張 係 数αと 温 度 変 化∆T の 積 で 決 定 さ れ る ．こ の 支 配 方 程 式 は ，マ ク ロ 一 般 化 ひ ず みE˜ を データとして，未知の擾乱変位u∗について面内周期性 の 拘 束 条 件 下 で 解 く べ き 問 題 と なって い る ．な お ，こ こでは簡単のために等方的な熱膨張変形 を仮定してい る が ，異 方 性 を 想 定 し て も 一 般 性 は 失 わ れ な い ．
一 方 ，マ ク ロ 一 般 化 応 力M˜ は ，次 式 で 定 義 さ れ る ． M˜ =
∫ h/2
−h/2
zT|z3→y3
( 1
l1l2
∫ l1/2
−l1/2
∫ l2/2
−l2/2
σdy1dy2
) dy3
=
∫
zTΣ(z3)dz3 (9)
上 式 のΣは マ ク ロ 応 力 で あ り，次 式 の よ う な 面 内 方 向 の 平 均 化 に よ り 定 義 さ れ る ．
Σ(z3)= 1 A
∫ l1/2
−l1/2
∫ l2/2
−l2/2σ(y1, y2, z3)dy1dy2 (10) こ こ で ，Aは 面 内 方 向 の 断 面 積A=l1l2で あ り，Vは 面 内ユニットセルの体積である．また，ここで定義したz
は 上 記 で 定 義 し た˜zと 異 な り，面 内 平 均 を と る 仮 定 で
−z1/2,−z2/2の 項 が0と なって い る[12, 13]．す な わ ち ，
z=
1 0 0 z3 0 0 0 0
0 1 0 0 z3 0 0 0
0 0 0 0 0 0 0 0
0 0 1 0 0 z3 0 0
0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 1
(11)
のよ うに定義されている．式(8)の第2，3式を式(9)に 代 入 す る と ，均 質 化 板 剛 性 をD˜ と し て ，次 式 の よ う な マ ク ロ 一 般 化 応 力M˜ と マ ク ロ 一 般 化 ひ ず みE˜ の 関 係 式 を 導 出 す る こ と が で き る ．
M˜ =D˜ (
E˜−∆T ˜Eth)
(12) こ こ で ，E˜thは 単 位 温 度 変 化 あ た り に 生 じ る マ ク ロ 一 般 化 熱 ひ ず み で あ り，面 内 ユ ニット セ ル の ト ポ ロ ジ ー に 応 じ て 決 定 さ れ る も の で あ る ．均 質 化 板 剛 性 と マ ク ロ 一 般 化 熱 ひ ず み は ，等 方 性 材 料 の 仮 定 を 用 い て 扱 わ れ る こ と が 多 い が ，本 研 究 で は 次 式 の よ う に 定 義 す る こ と で ，完 全 異 方 的 な 挙 動 を 許 容 す る も の と す る ．
D˜ =
A B K
BT D L KT LT H
=
A11 A12 A31 B11 B12 B13 K11 K12
A12 A22 A23 B21 B22 B23 K21 K22
A31 A23 A33 B31 B32 B33 K31 K32
B11 B21 B31 D11 D12 D31 L11 L12
B12 B22 B32 D12 D22 D23 L21 L22
B13 B23 B33 D31 D23 D33 L31 L32
K11 K21 K31 L11 L21 L31 H11 H12
K12 K22 K32 L12 L22 L32 H12 H22
(13)
E˜th={
E˜1th E˜2th E˜th3 E˜th4 E˜th5 E˜th6 E˜th7 E˜8th}T
(14) ま た ，本 研 究 で は 熱 伝 導 解 析 は 行 わ ず 常 に 一 様 温 度 を 与 え る こ と と し ，式(8)の 第4式 で 用 い た 温 度 変 化∆T を そ の ま ま マ ク ロ 構 造 に も 適 用 す る ．な お ，式(12)は マ ク ロ 一 般 化 熱 応 力M˜thを 用 い て 次 式 の よ う に 表 す こ と が で き る ．
M˜ =D ˜˜E−∆T ˜Mth (15) こ こ で ，マ ク ロ 一 般 化 熱 応 力 は ，M˜th=D ˜˜Ethと 定 義 し た ．な お こ の 構 成 関 係 は ，マ ク ロ 材 料 特 性 を 定 義 す る こ と な く，マ ク ロ 変 形 特 性 の 異 方 性 や 面 内 お よ び 面 外 変 形 の 連 成 作 用 を 表 現 し て い る こ と に 注 意 さ れ た い ．
本研究では，面内ユニットセルに数値平板試験[12,13]
を 行 う こ と に よ り，均 質 化 板 剛 性 D˜，並 び に マ ク ロ 一 般 化 熱 応 力M˜thを 算 出 す る ．詳 し く は 付 録A，あ る い は 文 献[12, 13]を 参 照 さ れ た い ．．
x
1
2
3
0
e1 0
e2 0
e3
1
3
2 e1
e2
e3
C0
C
y
z
Global
coordinate system Initial
configuration Initial
coordinate system
Co-rotational
coordinate system Co-rotational configuration
Deformed configuration
Fig. 4 Configurations and coordinate systems in co- rotational formulation
3. 共回転理論に基づくマクロ構造解析 手法
本 章 で は ，マ ク ロ 構 造 に 対 し て 大 変 位・大 回 転 問 題 を 取 り 扱 う た め の ，Felippaら[18]が 提 案 し た 共 回 転 理 論 と ，第2章 で 導 出 し た 均 質 化 板 剛 性 と 一 般 化 熱 応 力 の 適 用 法 に つ い て 解 説 す る ．共 回 転 理 論 は 第1章 で 触 れ た よ う に ，物 体 の 全 運 動 を 剛 体 運 動 と 弾 性 変 形 と 分 け て 考 え ，弾 性 変 形 に よ る 変 位 ，回 転 量 を 用 い て ひ ず み，応力を計算し，つり合い条件を満たす手法である．
そ の た め ，各 物 理 量 に つ い て ，全 運 動 に よ る 値 と ，剛 体 運 動 に の み ，あ る い は 弾 性 変 形 に の み 依 存 す る 値 と の関係を把握することが必要となる．そこで以下では，
まずFelippaら[18]の方法に倣い，全運動を剛体運動と
弾 性 変 形 に 分 解 す る た め 配 置 と 座 標 系 を 設 定 し ，物 理 量 同 士 の 関 係 を 示 す．次 に ，均 質 化 板 剛 性 を 用 い て 内 力 ベ ク ト ル を 導 出 す る 方 法 に つ い て 述 べ る ．な お ，本 研 究 の 定 式 化 で は 直 交 座 標 系 を 採 用 す る ．
3.1 配 置 と 座 標 系 の 設 定 共 回 転 理 論 で は 変 形 を 記 述 す る 際 に ，初 期 配 置（initial configuration），共 回 転 配 置（co-rotaitonal configuration），現 在 配 置（deformed configuration）の3つの配置を定義する（Fig.4）．初期配 置 は 運 動 す る 前 の 配 置 で あ り，現 在 配 置 は 運 動 後 の 配 置である．共回転配置は，これら二つの中間に位置し，
全運動のうち 剛体運動のみ生じさ せた仮想的な配置で あ る ．共 回 転 配 置 か ら 現 在 配 置 ま で の 変 動 が 弾 性 変 形 に よ り 生 じ た も の と 認 識 す る こ と が で き る ．ま た ，全 体 座 標 系 の ほ か に ，初 期 配 置 ，共 回 転 配 置 に つ い て そ れぞれ局所座標系を設ける（Fig.4）．これらの局所座標 系 を 設 定 し ，共 回 転 配 置 と 現 在 配 置 を 同 じ 局 所 座 標 系 で 議 論 す る こ と で ，物 体 の 全 運 動 か ら 剛 体 運 動 に よ る 影 響 を 取 り 除 く こ と が で き る ．
各 座 標 系 の 設 定 の 仕 方 は さ ま ざ ま で あ り，要 素 形 状 に よって も 変 化 す る が[22]，こ こ で は 本 研 究 で 採 用 し た3節 点3角 形 要 素 に 限 定 し 議 論 を 進 め る こ と に す る. 初 期 配 置 に お け る 局 所 座 標 系（ 以 降 ，初 期 座 標 系 ）の 原 点 をC0，ま た そ の 座 標 値 をx0C0と し ，各 初 期 座 標 軸
の 基 底 ベ ク ト ル をe01，e02，e03 と す る ．ま ず，原 点C0を 各有限要素の重心に設定し，一つ目の基底e01を有限要 素 の 任 意 の1辺 と 平 行 に な る よ う に 設 定 す る ．次 に 要 素の面に対して垂直になるようe03を設定し，最後のe02 に つ い て は 既 に 定 義 し た2軸 と 垂 直 に 交 わ る よ う 定 義 す れ ば ，す べ て の 座 標 軸 が 得 ら れ る ．こ れ ら を 数 式 で 表 す と 次 式 の よ う に な る
x0C0=1 3
∑3
a=1
x0a (16)
e01= x021
∥x021∥, e03= x021×x031
∥x021×x031∥, e02=e03×e01 (17) こ こ で ，x021=x02−x01で あ り，x0a，(a=1，2，3)は ，各 節 点 の 全 体 座 標 系 に お け る 座 標 値 で あ り，aは 各 要 素 ご と に 設 定 さ れ る 要 素 番 号 で あ る（Fig.4中 の1, 2, 3に 該 当 ）．ま た ，全 体 座 標 系 か ら 初 期 座 標 系 へ の 回 転 行 列 T0は
T0=[e01e02e03]T (18) のようになる．共回転配置における局所座標系（以後，
共 回 転 座 標 系 ）に つ い て は ，変 形 後 の 各 節 点 座 標 値 xa=x0a+ua，(a=1，2，3)を用いて，同様の手順で原点C と そ の 座 標 値xC，各 座 標 軸e1，e2，e3な ら び に ，回 転 行 列Tを 設 定 す る ．
最後に本章で用いる各変数の表記法についてまとめ て お く．各 座 標 系 に 関 す る 変 数（ 例 と し てxを 取 り 上 げ る ）に つ い て ，xの よ う に 上 線 を つ け た も の は 局 所 座 標 系 で の 値 で あ り，上 線 の つ い て い な い も の は 全 体 座 標 系 の 値 で あ る ．ま た ，変 数 に つ け る 上 付 き 文 字 に は，0，R，eの3つがあり，x0は変形前の変数を，xRは 剛 体 運 動（ 局 所 座 標 系 の 回 転 ）の み の 影 響 を 考 慮 し た 変 数 を ，xeは 弾 性 変 形 の み の 影 響 を 考 慮 し た 変 数 を 表 している．これらに対して特に上付き文字のないxは，
全 運 動 後 の 変 数 を 表 し て い る ．こ こ で ，局 所 座 標 系 に は2種類設定したが，x0は初期座標系，x，xR，xeは共 回 転 座 標 系 に よ り 表 し た 変 数 で あ る ．さ ら に ，本 章 で は マ ク ロ ス ケ ー ル に 絞って 議 論 を 展 開 す る た め ，2章 で 定 義 し た ミ ク ロ 座 標 系y，zは 用 い な い ．な お ，変 位 場 に つ い て は 第2章 で 用 い た ミ ク ロ 変 位 場uと 同 じ 表 記 を 用 い て い る が ，こ こ で のuは マ ク ロ 変 位 場 で あ る こ と に 注 意 さ れ た い ．
3.2 各 配 置 の 自 由 度 と 各 変 数 の 変 分 本 節 で は ， 各配置にて定義される節点変位と回転について説明し，
それ らの変分をとることで各変数 の微小増分量の関係 式 を 示 す．具 体 的 に は ，全 体 座 標 系 で 定 義 さ れ る 全 節 点 変 位 ，回 転 と ，共 回 転 座 標 系 で 定 義 さ れ る 弾 性 変 形 の み に よ り 生 じ る 全 節 点 変 位 ，回 転 と の 関 係 を 示 す こ と を 目 標 と す る ．ま ず，全 体 座 標 系 に お け る 弾 性 変 形 に 起 因 す る 変 位ueの 定 義 を 記 す．
uea=xa−xRa =( ua+x0a)
−(
x0C0+uC+xRaC) (19) また，変位ueaを図示したものをFig.5に示す．ここで，uC
x
a
a
C0 C
y
z
uC
ua
0
xa 0 0
xC
R
xaC e
ua
Initial configuration
Co-rotational configuration
Deformed configuration Global
coordinate system
Fig. 5 Illustration of elastic displacement ue
はC0からCまでの変位を表しており，またxRaC=xRa−xC
である．上式は全運動後 の節点座 標xaから，剛体 運動 の影響のみを受ける節点座標xRa を差し引くことで，弾 性 変 形 に よ る 変 位 を 導 出 し て い る ．共 回 転 座 標 系 に お け る 弾 性 変 形 変 位ueaは ，次 式 の よ う に な る ．
uea=Tuea=T(
ua+x0a−x0C0−uC
)−x0aC0 (20)
こ こ で ，上 式 の 右 辺 第2項 は ，局 所 座 標 系 に お け る 座 標 値 は 剛 体 回 転 に よって 変 化 し な い こ と を 表 す 次 式 を 用 い て 導 出 し た ．
x0aC0=xRaC0 (21) 次 に 回 転 自 由 度 に つ い て の 関 係 を 示 す．ま ず，各 節 点aご と に 定 義 さ れ る 全 回 転 行 列Raと 共 回 転 座 標 系 に お け る 弾 性 変 形 に 起 因 す る 回 転 行 列Rea と の 関 係 を 次 式 の よ う に 表 す．
Rea=RaR0T (22) 上 式 は ，全 回 転 行 列Raか ら 剛 体 回 転 に よ る 回 転 行 列 R0 の 逆 行 列 を 乗 ず る こ と で 弾 性 変 形 に よ る 回 転 行 列 を 導 出 し て い る ．共 回 転 座 標 系 で は 次 式 の よ う に 表 さ れ る ．
Rea=TReaTT=TRaT0T (23) こ こ の 式 展 開 で は 式(22)と ，3.1節 で 定 義 し たT0とT の 関 係 を 表 す 次 式 を 用 い て 導 出 し た ．
T=T0R0T (24) 一 方 ，回 転 自 由 度 を ベ ク ト ル 表 記 に す る た め に ，回 転 行 列Rに つ い て 以 下 の 擬 ベ ク ト ル を 用 い る ．
θ= θ
2 sinθN (25)
N={R(3,2)−R(2,3)R(1,3)−R(3,1)R(2,1)−R(1,2)} (26)
θ=∥N∥ (27)
上 記 の(3,2)な ど の 数 字 は 行 列 の 成 分 を 表 し て い る 。 こ れ をReaに も 適 用 さ せ る こ と で ，弾 性 変 形 に よ る 回 転 行 列 の 擬 ベ ク ト ルθeaを 導 出 す る こ と が で き る ．
本 研 究 で は ，3次 元 の 有 限 回 転 問 題 も 解 析 対 象 と な る た め ，回 転 自 由 度 の 更 新 の 際 に は 回 転 擬 ベ ク ト ル を 加 算 的 に 用 い る こ と は で き な い ．そ こ で ，3次 元 微 小
増 分 回 転 ベ ク ト ルδωeaを 用 い る こ と で 次 式 の よ う に 回 転 行 列 を 更 新 す る[17]．
Ra，new=R(δωea)Ra，old (28) R(δωea)=I+sin(δωea)
δωea
spin(δωea) +1
2
[sin(δωea/2) δωea
spin(δωea)2 ]
(29)
δωea=∥δωea∥ (30)
ここで，この3次元微小増分回転ベクトルδωeaは，回転 行列の擬ベクトルの増分δθと次式の関係にある[17,21]．
δθea=Ha(θea)δωea (31) Ha(θea)=I−1
2spin(θea)+ηspin(θea)2 (32) こ こ で ，上 式 のηは 次 式 の よ う に 求 め る ．
η=1−12θeacot12θea
θea
2 ≃ 1
12+ 1 720θea
2 (33)
θea=θea (34)
ま た ，spin(a)は ，あ る ベ ク ト ルa={a1,a2，a3}に 対 し て 次 の よ う な 成 分 を 有 す る 行 列 で あ る ．
spin(a)=
0 −a3 a2
a3 0 −a1
−a2 a1 0
(35)
最 後 に ，共 回 転 座 標 系 に お け る 弾 性 変 形 に 関 す る 変 位・回 転 自 由 度 を 含 む 一 般 化 変 位 ベ ク ト ル 微 小 増 分 δpe={δue, δωe}と 全 体 座 標 系 の 全 一 般 化 変 位 ベ ク ト ル 微 小 増 分δd={δu, δθ}は ，
δpea=Λδd (36)
の よ う に 表 す こ と が で き る ．こ の 行 列Λに つ い て は 付 録B.1を 参 照 さ れ た い ．
3.3 均質化板剛性を用いた内力ベクトルの設定 本 研 究 で 対 象 と す る 材 料 に は ，局 所 座 標 系 に お い て 微 小 ひ ず み 理 論 の 枠 組 み で 線 形 弾 性 体 を 仮 定 す る ．こ の と き，3.1節で各要素ごとに設定した共回転座標系を面内 ユ ニット セ ル の 支 配 方 程 式 の 記 述 に 用 い る と ，一 般 に 面内周期性の定義との不整合を生じる（Fig.6）．そのた め ，Fig.6に 示 す よ う に ，運 動 前 の マ ク ロ 構 造 要 素 に つ い て 初 期 座 標 系 の ほ か に 面 内 周 期 性 と 整 合 し ，均 質 化 板 剛 性 を 適 用 で き る よ う な 座 標 系（ 初 期 ユ ニット セ ル 座 標 系 ）を 設 定 す る ．こ の 初 期 ユ ニット セ ル 座 標 系 の 座 標 軸˜e01，˜e02，˜e03，お よ び そ の 回 転 行 列T˜0=[˜e01˜e02 ˜e03]Tは ， マク ロ構造の初期形状を平板として扱う場合には全体 座 標 系 の 原 点 をC0に 水 平 移 動 さ せ た 座 標 軸 を 採 用 す る こ と で 上 記 条 件 を 満 た す こ と が で き る ．こ の 初 期 ユ ニット セ ル 座 標 系 を 剛 体 回 転 さ せ た 座 標 系（ ユ ニット セ ル 座 標 系 ）へ の 回 転 行 列T˜ =[˜e1 ˜e2 ˜e3]Tは 式(24)と 同 様 に 次 式 で 定 義 さ れ る ．
T˜ =T˜0R0T (37)
Initial coordinate system
(example for NOT fulfilling in-plane periodicity)
Initial unit cell coordinate system (example for fulfilling in-plane periodicity) In-plane unit cell
Macrostructure y1
y2
y3
Macrostructure
Fig. 6 Setting for unit cell coordinate system
ま た ，共 回 転 座 標 系 に お け る 独 立 変 数xと ，ユ ニット セ ル 座 標 系 に お け る 独 立 変 数˜xの 関 係 は ，変 位 に 依 存 し な い 回 転 行 列R˜0を 用 い て 次 式 の よ う に 表 す こ と が で き る ．
˜x=R˜0x (38)
R˜0=T˜0T0T (39) 次 に ，ユ ニット セ ル 座 標 系 に て 内 力 ベ ク ト ル は 次 式 の よ う に 計 算 さ れ る ．
f˜=˜k ˜pe−f˜th (40) この式中の要素剛性行列˜kは均質化板剛性D˜ を用いて，
˜k=
∫
B˜TD ˜˜BdΩele (41) と表される．ここで，B˜ は，ユニットセル座標系におけ る 弾 性 変 形 に 起 因 す る 変 位・回 転 自 由 度p˜eと マ ク ロ 一 般 化 ひ ず みE˜ を 結 び つ け る 変 換 行 列 で あ り，E˜ =B ˜˜pe と な る ．ま た ，Ωeleは 各 マ ク ロ 構 造 要 素 の 全 領 域 を 指 している．加えて式(40)のf˜thは，マクロ一般化熱応力 を 用 い て 次 式 の よ う に 定 義 し た ．
f˜th=
∫
B˜TM˜thdΩele (42) さ ら に ，式(40)は 次 式 の よ う に 共 回 転 座 標 系 に お け る 内 力 ベ ク ト ル に 書 き 換 え る こ と が で き る ．
f=R˜(6)Tf˜=R˜(6)T˜k ˜R(6)pe−R˜(6)Tf˜th
=k pe−fth (43) こ こ で ，pe={ue,θe}は 共 回 転 座 標 系 に お け る 弾 性 変 形 に 起 因 す る 変 位・回 転 自 由 度 で あ り，前 節 で 導 出 し た ． また，R˜(6)=diag[
R˜0R˜0 R˜0 R˜0R˜0R˜0
]と定義している．最 後 に ，kは 次 式 の よ う に 定 義 し た ．
k=
∫
R˜(6)TB˜TD ˜˜B ˜R(6)dΩele=
∫
BTDBd˜ Ωele (44) こ こ で ，B=B ˜˜R(6)で あ る ．
上記の共回転座標系での内力ベクトルを全体座標系 で 表 す た め に ，次 の 関 係 式 を 用 い る ．
δdTf=δpeTf (45)
こ の 式 は ，全 体 座 標 系 の 仮 想 仕 事 と 共 回 転 座 標 系 の 仮 想 仕 事 が 等 し い こ と を 示 し て い る ．式(36)よ り，上 式 の 右 辺 を 展 開 す る と ，
δpeTf =δdTΛTf (46)
と な り，
f =ΛTf (47)
な る 関 係 式 が 導 か れ る ．
最 後 に ，こ の 内 力 ベ ク ト ル に つ い て 変 分 を と る と ， 次 式 の よ う に な る ．
δf =kTδd (48)
本 式 中 の 要 素 ご と に 定 義 さ れ る 接 線 剛 性 行 列kT に つ い て は 付 録B.2を 参 照 さ れ た い 。
4. マルチスケールトポロジー最適化
本 章 で は ，最 適 化 計 算 に お け る 設 計 変 数 の 設 定 と そ の 設 計 変 数 を 用 い て 表 さ れ る 面 内 ユ ニット セ ル の 有 効 弾 性 係 数 に つ い て 説 明 し た 後 ，変 位 量 制 御 問 題 に お け る 目 的 関 数 と 制 約 条 件 式 を 設 定 す る ．ま た ，本 研 究 で は感 度解析に基づく最適化アルゴリズムを適用するた め ，各 目 的 関 数 に つ い て の 解 析 的 感 度 を 導 出 し ，そ の 検 証 を 行 う．
4.1 設計変数および有効弾性係数の設定 設計領 域 と な る 面 内 ユ ニット セ ル は2種 数 の 線 形 弾 性 材 料 か ら な る 複 合 構 造 と し ，そ の 有 限 要 素 モ デ ル の 各 要 素 に おけ るどちらか一つの材料の体積分率を設計変数とす る ．こ の 設 計 変 数 を0≤si ≤1の 範 囲 で 変 化 さ せ る こ と で ，面 内 ユ ニット セ ル の 材 料 配 置 を 決 定 す る こ と に す る ．
熱 変 形 を 考 慮 す る 場 合 は ，ヤ ン グ 率 だ け で な く 熱 膨 張 係 数 も 設 計 変 数 に 依 存 す る ．本 研 究 で は ，熱 膨 張 係 数αと ヤ ン グ 率Eを 乗 じ た 熱 応 力 係 数β(=αE) [23]を 導入し，この係数について設計変数を用いて内挿する．
ま ず 式(8)第2式 で 与 え た ミ ク ロ 応 力 の 式 に つ い て 次 式 の よ う に 展 開 す る ．
σ=C( ε−εth)
=Cε−Cεth
=ECε−βC∆Tψ (49) こ こ で ，Cは 弾 性 係 数 行 列 の ヤ ン グ 率 に 依 存 し な い 項 である．この熱応力係数は，設計変数siを用いてRAMP 法[23, 24]に よ り 次 式 の よ う に 内 挿 す る 。
E(si)= (
1− si
1+qE(1−si) )
E1+
( si
1+qE(1−si) )
E2 (50) β(si)=
(
1− si
1+qβ(1−si) )
β1+
( si
1+qβ(1−si) )
β2 (51) 詳 細 に つ い て は 文 献[23]を 参 照 さ れ た い ．
4.2 変位量制御問題の設定 温度変化が生じたと き の マ ク ロ の 変 形 を 制 御 す る た め に ，指 定 し た 節 点 の 変 位 量 最 大 化 問 題 を 考 え る ．指 定 す る 節 点 番 号 をkと
し た と き の ，目 的 関 数 ，制 約 条 件 を 以 下 の よ う に 設 定 す る ．
min f (si)=fkUk (52)
subject to h(si)= (∫
V
sidV−V0
)2
−δ0<0 (53) ここで，Ukは指定した節点自由度の箇所だけ値を持つ 仮 想 荷 重 ベ ク ト ルFvと 各 要 素 ご と に 設 定 さ れ る 非 ゼ ロ の 仮 想 荷 重 係 数fkvを 用 い て 以 下 の よ う に 表 す．
fkvUk=FvTU (54)
Fv=( f1v,f2v,· · ·,fkv,fk+1v ,· · ·)T=(0,0,· · ·,fkv,0,· · ·)T (55) U=(U1,U2,· · ·,Uk,Uk+1,· · ·)T (56) 上記の仮想荷重ベクトルFvと節点変位ベクトルUの内 積 を と る こ と で ，任 意 節 点 の 変 位 量 の み を 抽 出 す る こ と が で き る ．ま た ，指 定 す る 節 点 は 複 数 選 ぶ こ と も で き ，そ の 場 合 に は 節 点 に 応 じ て 荷 重 係 数 の 絶 対 値 の 値 を 分 布 さ せ る こ と も 可 能 で あ る が ，本 解 析 で は fkv =1 に固定する．最後に，制約条件式(53)は，ある微小値δ0
を 用 い て ，体 積 一 定 条 件 と 等 価 な 不 等 式 制 約 条 件 式 を 与 え て い る ．こ こ で ，Vは 面 内 ユ ニット セ ル の 全 体 積 ， V0は 初 期 の 面 内 ユ ニット セ ル のphase-2の 体 積 で あ る ． 次に，目的関数の解析的感度を随伴法を用いて導出す るために，つり合い点の残差ベクトルr=Fext−Fint=0 と 随 伴 ベ ク ト ルλを 用 い て ，随 伴 関 数 を 次 式 の よ う に お く．
f∗(si)=FvTd−λTr (57) ここで，Fintは全体内力ベクトルであり，式(47)で定義 し たf を 全 マ ク ロ 構 造 要 素 に つ い て ア セ ン ブ リ す る こ と で 得 ら れ る ．ま た ，Fextは 全 体 外 力 ベ ク ト ル で あ る ．
式(57)の設 計変 数siにつ いて の導 関数 は次 式 の よ う に な る ．
∂f∗(si)
∂si
=dFvT dsi
U+FvT∂U
∂si
−λTdr dsi
=FvT∂U
∂si
−λT
(∂r
∂U
∂U
∂si
+ ∂r
∂si
)
=(
FvT−λTKT
)∂U
∂si −λT∂r
∂si (58)
ここで，KTは式(48)，(123)で定義した要素ごとの接線
剛性行列を全体系にアセンブリしたものである．また，
随 伴 ベ ク ト ルλは 任 意 で あ る た め ，下 記 の 関 係 を 満 た す よ う に 設 定 す る ．
Fv=KTλ (59)
そ し て ，最 終 的 に 式(58)は 次 式 の よ う に な る ．
∂f∗(si)
∂si
=−λT∂r
∂si
(60)
こ こ で ，∂r
∂si
の う ち ，設 計 変 数siに 陽 的 に 依 存 す る 項 は対応する要素の均 質化板剛性とマクロ 一般化熱応力
![Fig. 4 Configurations and coordinate systems in co- co-rotational formulation 3. 共回転理論に基づくマクロ構造解析 手法 本 章 で は ，マ ク ロ 構 造 に 対 し て 大 変 位・大 回 転 問 題 を 取 り 扱 う た め の ， Felippa ら [18] が 提 案 し た 共 回 転 理 論 と ，第 2 章 で 導 出 し た 均 質 化 板 剛 性 と 一 般 化 熱 応 力 の 適 用 法 に つ い](https://thumb-ap.123doks.com/thumbv2/123deta/7350986.2437008/4.892.463.810.108.299/FigConfigurations共回転理基づくマクロ構造解析手法本章マクロ構.webp)
![Table 2 Accuracy of sensitivities for homogenized plate sti ff ness D˜ α,β ∆ D˜ /∆ s i ∂ D˜ /∂ s i Error(%) A 1 , 1 [N / mm] 1](https://thumb-ap.123doks.com/thumbv2/123deta/7350986.2437008/9.892.88.496.166.482/table-accuracy-sensitivities-homogenized-plate-sti-ness-error.webp)
