1．は じ め に 自動車，建設機械，船舶，航空機等の大型機械から家 電，コンピュータ，マイクロデバイス等の小型機器に至 るあらゆる機械製品の開発において，環境配慮のための 低燃費化，コスト低減やユニバーサルデザイン，市場 ニーズへの迅速な対応等，その開発要求は非常に厳しい ものとなっている。そのため，機械の構造設計において は軽量化，合理化や特性の改良の短期間での達成が必須 となり，それらを効率的に実現するためのツールとして 有限要素法を中心とする数値シミュレーション技術が設 計の上流から下流，システム全体から細部構成部品に至 る様々なフェーズで幅広く利用され，設計の一プロセス を担うようになってきている。GUI (Graphical User Inter-face) の進歩も相まって，携わるエンジニアも一部の専門 技術者から，それを専門としない技術者，設計者へと裾 野が広がってきている。それに伴い，従来からの特性の 評価や予測のための構造解析に加え，感度解析や構造解 析の逆問題となる構造最適化解析も定着しつつあり，自 動設計も視野に入れたより効率的な構造決定が可能に なってきている。 構造最適化は設計変数で分類すると寸法と形状とトポ ロジーに，又は離散系（パラメトリック）と分布系（ノ ンパラメトリック）に，一方，目的関数で分類すると， スカラー関数（単目的）とベクトル関数（多目的）のよ うにも分類される。市販の汎用構造解析コードに組み込 まれているもののほとんどは種々の問題に対する汎用性 を考慮した離散系の解法である。一方，本稿で示す解法 は分布系の形状最適化に位置付けられ，多自由度な形状 決定が可能な点に特徴がある。連続体を対象にした寸 法，トポロジー最適化は幾何形状を固定して解析を行う が，形状最適化は幾何形状そのものを設計変数として扱 うため，広まりつつある三次元 CAD との直接リンクが 可能であるところに特徴があり，発展が期待されてい る。 本稿では構造最適化技術の中の形状最適化に関して， まず，その概要を紹介し，続いて，畔上により提案さ れ1) ，著者らがこれまで開発してきた数値解法（力法と 称す）について，構造設計において重要な強度，剛性， 振動の各問題へ適用する場合の定式化，及び解法を例題 や自動車構造部品の解析例を交えて示すこととしたい。
構造設計における分布系形状最適化問題の汎用的数値解法
下 田 昌 利*
A General Numerical Solution to Distributed Shape Optimization Problems in Structural Designs
Masatoshi SHIMODA*In this paper, a numerical shape optimization method of continua is presented for typical strength, rigidity and vibra-tion problems in structural designs. As the strength problems, the minimizavibra-tion problem of maximum stress and the shape determination problems that achieve a given desired stress distribution are formulated. The rigidity problems in-volve the minimization problem of external work and the shape determination problems that achieve a given desired dis-placement distribution. Also, the vibration problems involve maximization of eigen frequency with mode tracking. Each problem is formulated and sensitivity functions are derived using the Lagrangian multiplier method and the material de-rivative method. The traction method, which is a shape optimization method, is employed to find the optimal domain variation that reduces the objective functional. The proposed numerical analysis method makes it possible to design op-timal structures efficiently. Examples of computed results are presented to show the validity and practical utility of the proposed method.
Vol. 37, No. 1, 2003
* 機械工学科 助教授
2．形状最適化の研究の概要 2.1 その狙い 要求仕様を満足する構造を二次元もしくは三次元 CAD を用いて，形状として具現化することが構造設計の目的 である。ある構造を設計する際，一般に，設計者はまず 過去の設計を基に，候補となる設計案を幾つか考え，要 求される仕様と設計案に対する性能の比較を行い，その 中から最良のものを選択する。構造がわずかなパラメー タの組み合わせから構成される場合は比較的問題は少な いが，実際には多数のパラメータが存在し，その候補を 考案することさえ難しい場合がほとんどである。また， 最適な構造が候補の中に含まれていることは稀である。 構造物の幾何学的分布形状を設計変数と考える形状最適 設計の場合は，その候補は無限に存在する。時間的余裕 があれば，経験や勘，過去の試験や有限要素法等による 解析結果に基づき，選択された設計案に構造変更が繰り 返され，改良が施される。しかし，現実には時間的制約 と計算コストの制約を受け，その変更にも限界がある。 こうした状況を改善するためには，計算速度と少ない構 造解析へのアクセスで最適な設計を見出す最適化理論が 必要となる。近年の有限要素法と計算機の発達は最適設 計のための支援ツールとして不可欠なものになっている。 従って，最適設計を見出す実用的な理論の構築がキーで あり，その実現により，効率的で，しかも従来の限界を 超えるような形状設計が期待できる。 2.2 本研究（力法）の位置付け 最適設計に関する研究の歴史を振り返ると，1600 年代 に Galileo は平等強さの概念を導入し，はりの最適形状 を求めている。計算機を用いた近代の形状最適設計の研 究となると，1973 年，Zienkievictz and Campbell の論文に 端を発する2) 。以来，形状最適化の研究はパラメータ化 された設計問題に対して最適化理論を適用し，得られた 感度と数理計画法により有限次元空間で解探索を行う離 散系の解法が主流であった。そこでは，設計境界を表現 する設計変数の数が制限されるため，自由境界の設計は 困難という問題は解決されないままであった。併せて， 形状最適化の研究では境界の滑らかさと有限要素のリ メッシュの問題も問題として顕在化していた。これらを 同時に解決する方法として力法が 1994 年に登場した。 力法は楕円型境界値問題における形状最適化問題の数 値解法として提案された。ヒルベルト空間での勾配法を 利用し，分布系の最適化理論から導出される領域変動の 支配方程式を基に最適形状を求める方法であり，変分法 や最適制御を起源とする関数空間から解関数を求める分 布系の形状最適化問題の解法として位置付けられる。現 在主流となっている有限次元空間で解探索を行う離散系 の解法とは異なり，設計変数（形状）を分布関数として 扱う分布系の最適化問題として定式化するため，自由度 を制限することなく多自由度な自由曲面の形状設計が可 能であること，滑らかな境界形状が得られること，リ メッシュが実用上不要であることを主な特徴とする。ま た，これまでの分布系の形状最適化問題に関する主な研 究はフランス応用数学者らによる理論的な研究が中心で あった3)。著者らはこの分布系の形状最適化理論を構造 設計へ応用し，汎用性，実用性のある数値解法として構 築することを目的として，研究を進めてきた4)–7)。 3．領 域 変 動 力法による求解の手続きは，問題の定式化から始ま り，ラグランジュ乗数法と物質導関数を利用して感度関 数である形状勾配関数（境界法線方向のベクトル関数） を導出し，その結果を領域変動量の解析に利用する。 本章では分布系の形状最適化問題を定式化する前準備 として，速度法8)による領域変動の表現方法を簡単に紹 介する。 領域変動は初期領域 WRn, n2, 3 を定義域とした 1 対 1 写像，Ts(X): XW xWsを用いて表現することが できる。s は変動の履歴を表すことにする。この時，領 域 Q における変動拘束を仮定すると，領域の微小変動は 次のように表すことができる。 TsDs(X)Ts(X)DsVO(|Ds|) (1) ここで，速度場 V は Ts(X)の Euler 導関数を与えており， 次のような連続関数として定義することができる。 VCQ{VC1(W : Rn)|V0 in Q} (2) 記号 O (|Ds|)は lim |Ds|Æ0O (|Ds|)/|Ds|0 を意味する。また， 汎関数の導関数は次のように得られる。汎関数 J が分布 関数 fsの領域積分 (3) で与えられている場合，Euler 導関数 J˙は式(1)の関係を用 いて次式で与えられる。 J sdx s f W
Ú
(4) ここで，nnniViを表している。ベクトル n は外向き 単位法線ベクトルを表す。G は領域 W の境界を表す。ま た，本稿のテンソル表示では Einstein 総和規約と偏微分 表示(・),i∂(・)/∂xiを使用する。分布関数において，(・˙ ) は Euler 導関数（物質導関数），(・) は Lagrange 導関数 （形状導関数）を表す8) 。 汎関数 J が分布関数 fsの境界積分 (5) で与えられる場合，Euler 導関数 J˙は次式で与えられる。 (6) ここで，k は領域が 2 次元の場合は曲率，3 次元の場合 は平均曲率を表す。以下に示す各問題の一連の定式化で はこれらを公式として使用する。 4．強 度 問 題 強度問題として，疲労や静強度の設計で基本となる最 大応力最小化問題と，強度上のヒューズの設定等の応用 が可能である与えられた目標応力分布を実現する形状を 求める形状同定問題（応力規定問題と呼ぶ）を考える。 後に示す剛性や振動の問題においても同様の手続きが可 能であるため，最も扱い難い最大応力最小化問題につい て，定式化から感度関数の導出までの過程を示す。 4.1 最大応力最小化問題 最大応力最小化問題において，強度の指標としては次 式で定義されるミーゼス応力 sMを取り上げ，体積制約 条件の下でその最大値を最小化する問題を定式化し，形 状最適化で重要となる感度関数を導出する。 (7) 4.1.1 目的汎関数の変換 最大応力を最小化するミニマックス問題の目的関数は 一般に次のように表される。 , (8) ここで，saは規準化のための定数を表す。 ミニマックス問題に潜む最大応力の局所性による微分 不可能と特異性の問題を回避するため，及び力法の適用 を考慮して分布系の最適化問題として定式化するために は目的汎関数を積分形式で表す必要がある。これらの問 題を同時に解決するため，次のような積分形式の KS 関 数9) を用いる。 (9) KS関数は定数 r が十分大きな値をとる時，関数 sM(x) の最大値を抽出することができる。本槁では r を最大測 度抽出定数と呼ぶことにする。実際には r5 から 200 程 度を適用した。同様の目的で KS 関数の代わりに lpノル ムを用いることも可能である。 4.1.2 問題の定式化と形状勾配関数の導出 Fig. 1のように，初期領域 W，境界 G¤∂W の線形弾性 体が変動して（変動 V により）領域 Ws，境界 Gs¤ ∂Wsと なるものと仮定する。Gdesignは変動可能な設計境界を表 す。物体力 f（テンソル表示では fi, i=1, …, n）と表面力 P（テンソル表示では Pi）はそれぞれ Ws，G1に作用して いるものとする。 測度としてミーゼス応力を考えた場合，体積および各 荷重ケースでの状態方程式を制約条件とする最大応力最 小化問題は式(9)を用いて次のように表される。 Given W , f in W , P on G1, MoR (10) find Ws(or V ) (11) KS x Mx dx a s ( ( ))s ln exp ( ) r s s r 1 Ê ◊ Ë Á ˆ ¯ ˜ Ï Ì Ô Ó Ô ¸ ˝ Ô ˛ Ô
Ú
W max ( ) x M a s Œ È Î Í ˘ ˚ ˙ W s s x s s s s s s s s s s M xx yy yy zz zz xx xy yz xz 2 2 2 2 2 2 2 1 2 6 {( ) ( ) ( ) ( )} ˙ { ( ) } , J s s i in s nd sÚ
f¢f f k n G G J sd s f G GÚ
˙J sdx s nd s sÚ
f¢Ú
f n W G G(12)
subject to MM00 (13) a(v,w)l(w) for all wU (14) また，変分ひずみエネルギーを与える双一次形式 a(v,w)， 外力による変分ポテンシャルエネルギーを与える一次形 式 l(w)は式(11)，式(12)で定義される。v（テンソル表示で は vi），w（テンソル表示では wi）はそれぞれ変位と変分 変位を表す。U は変位の拘束条件を満たす適当に滑らか な関数空間を表す。M, M0はそれぞれ体積とその制約値 を，Rは正の実数の集合を表す。 (15) (16) なお，e（テンソル表示では eijkl）は弾性テンソルを表す。 この問題に対するラグランジュ汎関数 L は w(x), L を それぞれ状態方程式と体積制約に対するラグランジュ乗 数として次のように表される。これにより，制約付き最 適化問題が無制約最適化問題に変換される。 (17) 簡単のため表面力の作用する境界は法線方向には変動 しないこと(nlVl0 or G )，物体力は領域内で一定である こと(f0)を仮定すると，ラグランジュ汎関数の領域変 動に対する物質導関数 L˙ は，領域変動に対する速度場 V を用いて次式のように表される。 (18) (19) ここで，n は外向き単位法線ベクトルを表す。 ラグランジュ汎関数 L の v,w および L に関する最適性 条件は不等式制約条件に対する停留条件を考慮すると次 のように表される。
a (v,w)l (w) for all wU (20)
(21) L(MM0)0, MM00, L0. (22)(23)(24) ここで，式(20)は状態方程式と一致した v の支配方程 式，式(21)は随伴変数 w に関する支配方程式（随伴方程 式）である。また，式(22)から式(24)は不等式制約条件式 に対する Kuhn–Tucker 条件式を与え，体積制約に関する ラグランジュ乗数 L の支配方程式になっている。 これらの条件で決定された v, w および L を用いると すれば，ラグランジュ汎関数の導関数は次式で与えられ ることになる。 L˙lG(V ) (25) ここで，速度場 V の一次形式 lG(V )は次式で与えられる。 (26) on Gdesign¤G \Gfix (27) なお，G は形状勾配関数と呼ばれ，設計境界上で与えら れる感度関数を表す。 導出した v の支配方程式(20)と体積制約に関するラグ ランジュ乗数 L の決定式(22)–(24)は問題に依存しないが， 随伴方程式と形状勾配関数，もしくは速度場 V の一次 形式は問題に依存し，目的汎関数や制約条件の選び方に より異なったものになる。 G e v w f w n A ijkl k l i j i i M a , , exp 1 r s s r ◊ ◊ Ê Ë Á ˆ ¯ ˜ Ê Ë Á ˆ ¯ ˜ Ï Ì Ô Ó Ô ¸ ˝ Ô ˛ Ô◊ L lG G V di i s ( )V G G
Ú
a A dx A v dx U M a M ij ij M a M ij ij k k s s ( , ) exp exp ¢ Ê ◊ Ë Á ˆ ¯ ˜ ◊ ¢ ◊ Ê Ë Á ˆ ¯ ˜ ◊ ◊ ◊ ¢ ¢ ŒÚ
Ú
v w v 1 1 s s r ∂s ∂s s s s r ∂s ∂s ∂s ∂s W W for all A a dx M a s s s s r expÊ ◊ Ë Á ˆ ¯ ˜Ú
W ◊n V dl l G L ¢(M M 0), VŒCQ ˙ ( ) ( , ) exp ( , ) , , exp L l a A dx a e v w f w A M a M ij ij ijkl k l i j i i M a s s ¢ ¢ ◊ Ê Ë Á ˆ ¯ ˜ ◊ ¢ ¢ ◊ ◊ ◊ Ê Ë Á ˆ ¯ ˜ Ï Ì Ô Ó Ô ¸ ˝ Ô ˛ ÔÚ
Ú
w v w v w 1 1 s s r ∂s ∂s s r s s r W G L L dx l a M a s ( , , , ) ln exp ( ) ( ) ( , ) ( ) W L L W v w x w v w M M 1 0 r s s ◊r Ê Ë Á ˆ ¯ ˜ Ï Ì Ô Ó Ô ¸ ˝ Ô ˛ ÔÚ
l f w di i P w di i s ( )w W G G WÚ
1Ú
a e v w dijkl k l i j s ( , )v w , , W WÚ
that minimize 1 r s s r ln exp M( ) a s dx x ◊ Ê Ë Á ˆ ¯ ˜ Ï Ì Ô Ó Ô ¸ ˝ Ô ˛ ÔÚ
W4.2 応力規定問題 与えられた目標応力分布を実現する応力規定問題では 任意の目標応力分布に対して，解が得られる保証はない ため，目的汎関数として，逆問題の定式化で一般に用い られる目標分布からの二乗誤差，すなわち，実際のミー ゼス応力分布 sM(x)と目標応力分布 ˜sM(x)を用いて，最小 化する目的汎関数を次式のように定義する。 (28) ここで，WC（変動後 WCs）は目標応力分布が与えられ る規定領域を表す。 最大応力最小化問題の場合と同様の手続きに従い，最 適性の条件を導出すると，式(21)に対応する問題依存の 随伴方程式は次式のようになる。
for all vU (29) また，同様に問題に依存する式(26),(27)に対応する速度 場 V の一次形式 lG(V )，及び形状勾配関数 G 1, G2 は式 (30)から式(32)のようになる。 (30) (31) (32) 5．剛 性 問 題 剛性問題としては，外力仕事最小化問題と変位規定問 題を取り上げることとする。導出の過程は省略し，最終 的に得られる問題固有の随伴方程式と形状勾配関数につ いて，その結果を示す。 5.1 外力仕事最小化問題 剛性の指標として，変位 v を用いて，目的汎関数に式 (16)で定義される外力仕事 l(v)を選んだ場合，外力一定を 仮定すると，l(v)の最小化は負荷点変位の最小化（剛性 の最大化）となる。強度問題の場合と同様の手続きによ り得られる随伴方程式と速度場 V の一次形式 lG(V )，形 状勾配関数 G はそれぞれ式(33)，式(34) (35)のように導出 される。 (33) (34) (35) ここで注目されるのは，目的汎関数を外力仕事に選ん だ場合，式(33)の随伴方程式と式(20)の状態方程式が一致 し，vw の関係が成立することである。このような問 題は自己随伴問題と呼ばれ，随伴方程式の解析が不要と なるため感度関数の計算が簡単になる。また，この問題 の場合，形状勾配関数 G はポテンシャルエネルギー密度 の 2 倍の値と等価となる。 5.2 変位規定問題 機構を考慮した構造設計への応用が期待できる与えら れた目標変位分布を実現する形状を求める問題（変位規 定問題と呼ぶ）を考える。4.2 の応力規定問題と同様， 目的汎関数に次式のような二乗誤差を定義する。 (36) ここで，GDは目標変位分布が与えられる規定境界，v˜iは 目標変位ベクトルの成分を表す。 この問題に対する随伴方程式は式(37)のように導出さ れる。 (37) また，変位規定境界および表面力の作用する境界が法 線方向に変動しないことを仮定すると，速度場 V の一次 形式場 lG(V )と形状勾配関数 G は次のように得られる。 (38) (39) G( e v wijkl k l, i j,f wi iL)n lG G V di i s ( )V G G
Ú
a v x v x v di i i Cs ( , )v w¢ 2Ú
{ ( ) ˜ ( )} ¢ W G {( ( ) ˜ ( )}{( ( ) ˜ ( )}v xi v xi v xi v x di D G GÚ
G(2 f vi ie v vijkl k l i j, ,L) n on Gdesign∫G G\ fix
lG G V di i s ( )V G G
Ú
a( , )v w¢ l( ) v¢ for all v¢ ŒU G2( )n , , e v wijkl k l i j f wi i L G1s s 2n M( ) ˜ ( )x Mx{
}
lG G V di i G V di i s cs ( )V 1 2 G G G WÚ
Ú
∂ a x x v v dx M M M ij ij k k cs ( , )v w¢ 2Ú
{
s ( ) ˜ ( )s}
◊∂s ◊ ◊ ¢ ∂s ∂s ∂ W sMx sMx dx Cs ( ) ˜ ( ){
}
Ú
2 W6．固有振動問題 固有振動数のシフトは機械の振動対策の基本である。 ここでは，体積制約を有し，初期に指定した r 次モード v(r)を追跡しながら，その振動固有値 l(r)を最大化する問 題について，同様の手順でその最適性の条件を示す。 この問題は前述の強度，剛性問題の定式化における目 的汎関数を特定の振動固有値l(r)に，状態方程式(14)を 次の固有振動方程式に置き換えることにより定式化され る。
a(v(r),w)l(r)b(v(r),w), for all wU (40)
ここで，双一次形式 b(v(r),w)は次式で定義される。 (41) なお，ここでの w は変分固有振動モード，r は密度を表 す。 強度， 剛性の場合と同様の導出手続きにより， La-grange 汎関数 L の v(r)，w，及び L に関する最適性条件 は次のようになる。
a(v(r),w)l(r)b(v(r),w), for all wU (42) a(v(r),w)l(r)b(v(r),w), for all v(r)U (43) b(v(r),w)1 (44) L(MM0)0, MM00, L0. (45)(46)(47) ここで，式(42)と式(44)は固有振動モード v(r)を決定するた めの固有振動方程式と正規化条件を与え，式(43) は随伴 固有振動モード w を決定することができる随伴固有振動 方程式を与えている。これらの式を比較すると，この問 題でも自己随伴関係 v(r)w が成立することがわかる。一 方，式(45)から式(47)は式(22)から式(24)と同様の不等式制 約条件式に対する Kuhn–Tucker 条件式である。 これらの条件で決定された v(r)，w と L を用いること にすれば，Lagrange 汎関数の導関数 L˙，または速度場 V の一次形式と形状勾配関数 G は次のように表すことがで きる。 (48) (49) 7．数 値 解 法 2章では数値解法である力法の特徴を紹介したが，こ こでは 4 章から 6 章で導出した形状勾配関数を用いた力 法による領域変動量の数値解析方法，及び汎用システム の概要を示す。 7.1 力法 力法はヒルベルト空間での勾配法を利用して目的汎関 数が減少するような領域変動量（速度場 V）を支配方程 式(50)に基づいて解く方法である。
a(V,w)l(w) for all wCQ (50) 式(50)の支配方程式は境界あるいは領域に負の形状勾 配関数G を外力として作用させた時の変位場として速 度場 V を解析することを示している。すなわち，力法で は領域変動が形状勾配関数を擬似弾性問題の外力として 作用させた時の変位場として求められる。従って，式(50) は通常の線形弾性問題の解法を用いて解くことができ， 一般に，有限要素法を用いる。 この支配方程式で決定された領域変動 V はラグラン ジュ汎関数 L を減少させることを確認してみる。状態方 程式と制約条件式が満たされている時，ラグランジュ汎 関数 L の摂動展開は次のように表すことができる。 DLlG(DsV )O(|Ds|) (51) ここで，式(50)を式(51)に代入して，弾性テンソル eijkl の正定値性に基づく a(v,w)の正定値性 (52) を考慮すると，Ds が十分小さい時，次の関係が得られ ることになる。 DLa(DsV,DsV )0 (53) この関係は式(50)で決定された速度場 V を用いて領域 を変形していけば，凸性が保証されている問題において， ラグランジュ関数 L は必ず減少する関係を与えている。 得られた領域変動量を用いて領域を更新する手続きを 繰り返すことにより，滑らかな収束形状が得られる。そ の数学的証明と理論の詳細については文献に譲る10) 。 7.2 形状最適化システム これまで力法を汎用化，実用化するため，汎用 FEM コードを利用したシステムを開発してきた。最適化のモ $a0: ( , )ax xa x2 for all xŒU G( e vijkl ( ) ,r k lv( ) ,r i jl r( )r v( )r iL)n lG G V di i s ( )V G G
Ú
b r v w dxi i s (v( ), )w r WÚ
ジュールを持たない一般的な FEM モジュールのみを有 する市販の構造解析コードにおいても，最適設計の解析 を可能にするものである。市販の汎用 FEM コードを一 種のサブルーチンとして状態方程式（応力解析と呼ぶ）， 随伴方程式（随伴解析と呼ぶ）および力法の支配方程式 (50)（速度解析と呼ぶ）を解いている。 随伴解析において外力として必要な法線方向表面力は 圧力荷重によって代替しているので，法線方向の計算は 不要である。応力解析と随伴解析，速度解析および形状 更新を繰り返すことにより，目的汎関数が最小化され， 最適形状を得ることができる。実行に際しては形状勾配 関数の相対性を利用し，1 回の領域変動量の大きさを調 整する領域変動係数 Ds を設定し，G に乗じることによ り 1 回の領域変動量を調整する。また，体積の制約値， 設計境界 Gdesignを設定する。なお，速度解析では応力解 析で設定される力学的境界条件とは異なる領域変動の拘 束条件（側面制約条件）を独立に設定する。システムの 流れを Fig. 2 に示す。 7.3 随伴変数の数値計算 形状勾配関数に含まれる状態方程式の解 v や v(r)は適 切な境界条件を与えることにより有限要素法により素直 に求めることができる。一方，随伴方程式の解 w を得る にはやや工夫を要する。ここではこれまでに示した各問 題に対する随伴方程式の具体的な数値解析方法を順に示 す。 まず，最大応力最小化問題における随伴方程式につい て，式(21)の随伴方程式をマトリックス表記すると式(25) のようになる。 (54) ここで，[B]Tはひずみ−変位マトリックスの転置，[D]は弾 性マトリックスを表す。領域内に の初期ひずみを外力として与えることにより随伴変位 w は求められる。また，初期ひずみを熱ひずみに変換して 与えることも可能である。 同様にして，応力規定問題の場合は 2{sM(x)s˜M(x)} · の初期ひずみを応力規定領域に，変位規定問題の 場合は 2{vi(x)-v˜i(x)}の表面力を変位規定境界に外力とし て与えることにより随伴変位 w を求めることができる。 7.4 制約条件の扱い これまで力法の実行において制約条件を考慮する方法 として，単一制約に対して有効な PID 制御の考えを用い る方法や複数制約にも対応できる感度マトリックスを利 用する方法を提示してきた。次章に示す計算結果では PID制御の考えを用いる方法11) を使用している。具体的 には，体積制約を満たすように決定されるラグランジュ 乗数 L は力G の中では一様な表面力とみなすことがで き，大きさ L の一様表面力を制御することにより体積制 約条件を満足させている。 8．計 算 結 果 本章では導出した強度，剛性，振動の各問題に対する ∂s ∂s M ij [ ] [ ][ ] { } [ ] [ ] exp B D B dx x B D A dx T T M a M ij s s W W
Ú
Ú
Ê ◊ Ë Á ˆ ¯ ˜ ◊ Ï Ì Ô Ó Ô ¸ ˝ Ô ˛ Ô 1 s s r ∂s ∂sFig. 2. Schematic flow chart of shape optimization system.
形状勾配関数を用いた解析例を示す。なお，いずれの問 題も初期の体積（面積）を制約値とした。また，二次元 問題の解析は平面応力状態を仮定し，数値解析には 4 節 点 1 次要素を，三次元問題には 8 節点 1 次要素を使用し ている。 8.1 強度問題 8.1.1 最大応力最小化（ソリッドアーム） 三次元問題への適用例として，圧縮荷重を受けるソ リッドアームの最大応力最小化問題の設定を Fig. 3 に示 す。応力解析の境界条件は Fig. 3(a)に示すように円孔周 辺を完全拘束し，他端に圧縮分布荷重 P を負荷した。速 度解析の制約条件は Fig. 3(b)に示すように円孔周辺およ び長手方向を滑り拘束し，幅方向と厚さ方法に領域変動 可能な制約を与えた。なお，数値解析は対称性を考慮し た二分の一モデルで行った。 Fig. 4に初期形状と得られた最適形状，および各ミー ゼス応力分布を示す。Fig. 5 に初期値で基準化した最適 化履歴を示す。目的汎関数と最大応力は最小化されてい ることが確認できる。最大応力は約 30% 減少し，円孔 の一部周辺を除きほぼ全域で均一応力分布になっている ことがわかる。 8.1.2 応力規定問題（形状同定問題） Fig. 6にノッチ平板の問題設定を示す。ノッチ近傍の 領域 WCに目標応力分布 s˜M(x)を与え，下面を設計境界 とした。与えた目標応力分布は Fig. 7 の形状について事 前に求めたものであり，この形状を正解（目標）形状と した。応力解析では Fig. 6(a)に示すように左端部を完全 拘束し，反対面に分布荷重を負荷した。速度解析では Fig. 3. Solid arm problem.
Fig. 4. Initial & optimal shape and stress distribu-tions (MPa).
Fig. 5. Iteration histories.
Fig. 6. Notched plate problem.
Fig. 6(b)に示すように応力規定領域，拘束，負荷境界お よび上辺は変動しないものとした。 Fig. 8に初期形状および得られた収束形状を示す。Fig. 9に応力規定領域（A から B）での応力分布を示す。Fig. 10に初期値で基準化した最適化履歴を示す。これらの結 果から面積は一定に保たれ，狙い通り目的汎関数はゼロ になり，応力分布も目標応力分布に一致していることが わかる。目標形状の Fig. 7 と比較するとやや波打ってい るが，目標とする応力分布を実現する形状が得られてい ることが確認できる。 8.2 剛性問題 8.2.1 外力仕事最小化問題（トーションアーム） Fig. 11に二次元トーションアームの最大変位最小化の 問題設定を示す。応力解析では Fig. 11(a)に示すように円 孔周辺を完全拘束し，他端に分布荷重を負荷した。速度 解析では Fig. 11(b)に示すように円孔周辺を固定し，両端 を滑り拘束した。 Fig. 12と Fig. 13 に初期形状と得られた最適形状，およ び 2 つの荷重ケースでの変形モードとミーゼス応力分布 を示す。Fig. 14 に初期値で基準化した最適化履歴を示 す。目的汎関数と最大変位は最小化されており，本解法 が有効に機能していることが確認できる。 比較のために，最大応力最小形状を求めた結果を Fig. 15に示す。剛性最大形状（外力仕事最小化）と強度最 大形状（最大応力最小化）とは必ずしも一致しないこと が確認できる。これは随伴問題のひずみ場の相違に起因 する。 8.2.2 変位規定問題（形状同定問題） 問題設定を Fig. 16 に示す。変位規定境界 GDを上面に 設定し，平面保持をするような一様な y 方向変位 v˜2を 与えた。応力解析および随伴解析では Fig. 16(a)に示すよ うに両下端部を完全拘束し，上面に分布荷重を負荷し た。速度解析では Fig. 16(b)に示すように変位モード規定 境界および荷重支持点は変動しないものとし，残りの全 境界を設計境界とした。 Fig. 17に初期形状と変形図を示す。なお，変形図は拡 大して示している。Fig. 18 に得られた形状と変形図を， Fig. 8. Shape identification result.
Fig. 9. Comparison of stress distributions along A–B.
Fig. 10. Iteration histories.
Fig. 19に初期値で基準化した最適化履歴を示す。これら の結果から目的通り変位規定境界の上面は一定に変位 し，目的汎関数もほぼゼロになっていることがわかる。 与えられた変位分布を実現する形状が得られ，変形モー ドのコントロールが達成されていることが確認できる。 Fig. 12. Initial shape.
Fig. 13. Optimal shape.
Fig. 14. Iteration histories.
8.3 固有振動問題 力法では自由曲面の設計が可能なため，鋳造や鍛造部 品の形状設計に対して特に有効である。ここでは自動車 部品への適用例を紹介する。 エンジンを車体に弾性支持するためのエンジンマウン トブラケットを対象に，固有振動問題の解析を行った。 初期形状と固有値解析の境界条件を Fig. 20(a)に示す。速 度解析 (Fig. 20(b)) ではこの円筒部分に変動拘束を与え， また，エンジンマウントブラケットの一部には平面を保 つ必要性から MPC (Multi Point Constraints) を用いて変動 制約を与えた。Fig. 21(a)に初期形状と振動モードを，(b) に体積一定条件で求めた基本固有振動数の最大化形状を 示す。この計算により，基本固有振動数が約 20% 増加 した収束形状が得られた。 9．お わ り に 本稿では構造設計で用いられている数値シミュレー ション技術の中の形状最適化において，開発してきた数 値解法（力法）を構造設計において重要で頻繁に遭遇す る強度，剛性，振動の代表的問題へ応用するための解法 について示した。併せて二次元，三次元の解析例を紹介 した。本方法により，強度，剛性，固有振動数の最大化 や応力分布，変位分布のコントロールを目的とする最適 形状設計（例えば，ヒューズ設計，機構設計）が可能に なる。 最適化手法は線形問題の範囲ではほぼ実際の構造へ適 用できるレベルになってきている。しかし，あくまで初 期構造の近傍で解構造を探索するものであるため，得ら Fig. 16. Table support problem.
Fig. 17. Initial shape and deformation mode.
Fig. 18. Shape identification result.
れる構造にも限界がある。基本となる初期構造の善し悪 しが結果に大きく影響することに常に留意する必要があ る。また，実際の構造設計問題を数理モデルに置き換え るため，いわゆる‘かたさ’を含んでいる。AI 等の技術 が応用され，‘やわらかさ’が考慮されてくると設計者 にとって使い勝手の良いツールになり，設計現場への普 及が加速されると思われる。また，実際の種々の制約条 件が多領域最適化問題として考慮されるにつれ，制約条 件や設計変数の数の増大による大規模最適化問題の効率 化や，並列分散処理等も要求され始めている。 非線形問題に関しても研究に進展が見られ始めている が，感度解析の計算負荷が大きいため，感度を近似的に 差分で求めてベーシスベクトル法と組み合わせる方法や 感度解析が不要な応答曲面法の研究が実用化が主流であ る。これらの方法は設計変数の数に実用上制限があり， 得られる最適構造は限定されたものならざるを得ない。 今後，効率的な感度解析の基礎理論と多自由度な設計変 数を扱える方法の研究の発展にも期待したい。 文 献 1) 畔上，機論，60-574, A(1994), 1479.
2) Zienkiewicz, O. C. and Campbell, J., in “Optimal Struc-tural Design”, Gallagher, R. H. and Zienkiewicz, O. C. (Eds.), John Wiley, New York (1973).
3) Cea, J., in “Optimization of Distributed Parameter Struc-tures”, Haug, E. J. and Cea, J., Vol. 2, Sijthoff & Noord-hoff, Alphen aan den Rijn (1981).
4) 下田，畔上，桜井，機論, 62-604, A (1996), 2831–2837 5) 下田，畔上，桜井，機論, 63-607, A (1997), 610–617. 6) 呉，畔上，下田，機論, 61-587, C (1995), 49–54.. 7) Ihara, H., Azegami, H. and Shimoda, M., in “Computer
Aided Optimum Design of Structures VI ”, WIT Press, Fig. 20. Engine mount bracket problem.
Southampton, (1999), 87–95.
8) Sokolowski, J. and Zolesio, J. P., in “Introduction to Shape Optimization, Shape Sensitivity Analysis”, Springer-Ver-lag, New York, (1991).
9) Kreisselmeier, G. and Steinhauser, R., in “Systematic Control Design by Optimizing a Vector Performance Index”, International Federation of Active Control
Sym-posium on Computer-Aided Design of Control Systems, Zurich, Switzerland, August 29–31, (1979), 113–117.
10) 畔上，海津，下田，片峯，第 1 回計算工学会講演
会論文集，Vol. 1, （1996 年 5 月).
11) Shimoda, M., Azegami, H., Sakurai, T., in “Proceedings of 9 th. International Conference on Vehicle Structural Mechanics and CAE”, Troy, April 4–6, SAE (1995), 223.
