勾配法と逐次二次計画法によるシェルの形状・板厚分布・位相最適化
Shape, Thickness Distribution, and Topology Optimization of a Continuum Shellusing the Gradient Method and Sequential Quadratic Programming◦ 渡辺 哲平 ∗1 ，藤田 慎之輔 ∗2Teppei WATANABE∗1 ，Shinnosuke FUJITA∗2∗1 北九州市立大学大学院 国際環境工学研究科 大学院生
Graduate Student, Faculty of Environmental Engineering, The University of Kitakyushu∗2 北九州市立大学大学院 国際環境工学研究科 准教授
Assoc. Prof., Faculty of Environmental Engineering, The University of Kitakyushu, Dr. Eng.
キーワード
キーワード：シグモイド関数; トポロジー最適化; シェル構造物
Keywords：Sigmoid function; topology optimization; shell structure.

1. 序論
シェルの形状や板厚分布は，力学的性能を大きく左右2. 設計問題の定式化するため，近年では，それらを設計変数とした最適化問

シェルの節点座標は (I+1) × (J+1) の制御ネットを有題を定式化し，構造最適化によって高い力学的性能を有する B スプライン曲面により離散化するものとする．

する形状や板厚分布を獲得する研究が行われている．形

加えて，各節点における板厚値も B スプライン曲面に状と板厚分布を直接的に設計変数として扱うと計算速度より離散化する．各制御点の x 座標を並べたベクトルをの低下や解の発散を招く恐れがあるため，パラメトリッ

q x ，y 座標を並べたベクトルを qy ，形状に関する制御点ク曲面により離散化する手法が提案されている 1) ．まの鉛直方向座標を並べたベクトルを qz ，板厚に関する制た，板厚分布最適化において節点板厚にシグモイド関数

御点の鉛直方向座標を並べたベクトルを qt とおく．

を乗じたトポロジー最適化手法も提案されている 2) ．既
往の研究では，トラス構造物の材料非線形解析や平面骨組みの幾何学的非線形解析に対して，剛性方程式を解く代わりに全ポテンシャルエネルギー最小化問題を解く数値実験が行われ，有効性が示されている 3–5) ．剛性方程
式を直接解く方法は，線形問題の場合には剛性方程式を1 度解くだけで解が得られるが，板厚が消失する場合においては，行列が特異になることが考えられる．一方，ポテンシャルエネルギーを最小化する方法は，線形方程()⊤q x = q x,00 · · · q x,0J · · · q x,I0 · · · q x,IJ()⊤qy = qy,00 · · · qy,0J · · · qy,I0 · · · qy,IJ()⊤qz = qz,00 · · · qz,0J · · · qz,I0 · · · qz,I J()⊤qt = qt,00 · · · qt,0J · · · qt,I0 · · · qt,I J

(1a)
(1b)
(1c)
(1d)

節点数を n，要素数を m とし，節点座標ベクトル r ならびに節点における板厚値 (節点板厚) を並べたベクトル tを次式で定義する.

式を直接解く代わりに，行列・ベクトル積の反復によって解が得られるため本研究の最適化問題に適した解法であると考えられる．本研究では，形状と板厚分布を B ス
プラインで離散化した上でシグモイド関数を板厚に乗じ()
r = (x y z)⊤
x = ( x1 · · · xn)⊤
y = (y1 · · · yn)⊤
z = z1 · · · zn()⊤
t = t1 · · · tn

(2a)
(2b)
(2c)
(2d)
(3)

ることで，連続最適化問題として定式化し，連続体シェル構造物に対して剛性方程式を解く代わりに全ポテン

このとき，節点座標並びに節点板厚は次式のように離散シャルエネルギー最小化問題を解き，形状と板厚分布，化される.

トポロジーの同時最適化を行う．さらに，全ポテンシャ
ルエネルギー最小化問題を複数の勾配法で解き，アルゴ

xk =I ∑
J∑q x,i j Bi,N (uks )B j,M (vks ),k = 1, · · · , n (4a)qy,i j Bi,N (uks )B j,M (vks ),k = 1, · · · , n (4b)i=0 j=0

リズムの違いによる解の収束性について考察する．

yk =I ∑J∑
i=0 j=0

日本建築学会情報システム技術委員会
第 44 回情報･システム･利用･技術シンポジウム論文集， 519-522 ， 2021 年 12 月，京都
Proceedings of the 44th Symposium on Computer Technology of Information ，Systems and Applications ， AIJ ， 519-522 ， Dec. ， 2021 ， Kyoto－ 519 －I ∑J∑zk =ここで，u は節点変位ベクトル， K は線形弾性行列であqz,i j Bi,N (uks )B j,M (vks ),k = 1, · · · , n(4c)i=0 j=0tk =I ∑J∑r∗ は支持点の節点座標ベクトル， r∗ini は r∗ の初期値でqt,i j Bi,N (uks )B j,M (vks ),k = 1, · · · , n(5)ある．tmin ，tmax は要素板厚の上下限値，qz,min ，qz,max ，qt,min ，qt,max は設計変数の上下限値である．

i=0 j=0(る．シェルの V は総体積，V ，Vmin は体積の上下限値，)vks は k 番目の節点の節点座標に対応したパ

ここで， 外力ベクトルを p とし，節点変位ベクトルを求ラメータ， Bi,N (u) は u 方向の N 次の B スプライン基底めることを考える．剛性方程式と全ポテンシャルエネル関数であることを表す．

ギーの最小化問題を，それぞれ式 (11) と式 (12) に示す．

ここで uks

板厚が非常に小さい値をとった場合，実際の施工を行Ku = pうことは困難になることを防ぐため，節点板厚に閾値t¯min を設け，t¯min を下回る際にはシェル形状に開口部をminimizeu開けることを考える．このとき，節点板厚は次のように

f (u) =(11)1 ⊤u Ku − p⊤ u2(12)

式 (11) を解くことと式 (12) を解くことは等価である．

表される．

tkt=p k
0iftk ≥ t¯miniftk < t¯min

式 (11) において tes = 0 で剛性マトリクスが特異になる

k = 1, · · · , n(6)ため式 (12) を用いて節点変位ベクトルを求める．

本研究で用いる有限要素解析は要素内で板厚は一定であ

3. 数値解析例

り，節点板厚は定義されない．そこで，節点板厚の平均3.1. 解析概要

値を次のように定義し，有限要素の要素板厚として用い球殻シェルを考える (図 1)．境界条件は図中に • で表しることを考える．

te =一辺が 30m の正方形平面を有するライズ 10m の裁断

p tei + p te j + p tek

3た四点を固定支持とする．形状並びに板厚分布はいずれ

e = 1, · · · , m(7)ei ，e j ，ek は要素 e の頂点における節点番号であり，p tei ，ijk
p te j ，p tek は e ，e ，e に対応する節点板厚である．しかし，式 (6) は te = t¯min で不連続となるため，感度解析を行うことができない．そこで，式 (6) にシグモイド関数を乗じたものを節点板厚として定義する．

も 3 次の B スプラインにより離散化し，制御点は 6 × 6格子状に 36 個配置する (I = J = 5)．節点数は n = 121，要素分割は m = 200 とする．外力は自重として 24kN/m3を作用させる．材料はコンクリートを想定する．初期のシェル厚は一様に 0.2m とする．図中の fini ，V として表した値は，初期形状におけるひずみエネルギーと体積を示す．本研究では，解析モデルに対して CG 法，

tetes =1 + exp{α(t¯min − te )}(8)

NewtonCG 法，加速勾配法 7) の 3 つの勾配法を式 (12)に適用し，収束性能の考察を行う．CG 法，NewtonCG

本研究では，設計変数 X に qz ならびに qt をとり，以下の最適化問題を解くことでシェルの板厚分布最適化，法については python の科学技術計算ライブラリである

scipy8) を用いる．

シェルの形状と板厚分布の同時最適化を行う．最適化計
算には pyOpt6) を用い，逐次二次計画法により解く．

minimizeXsubject toif
ifif

報告 H951 ⊤u Ku
2Vmin ≤ V ≤ Vr∗ = r∗initmin ≤ tes ≤ tmax e = 1, · · · , mqz,max ≤ qz,i j ≤ qz,max
i = 0, · · · , Iqt,max ≤ qt,i j ≤ qt,maxj = 0, · · · , J
f =X = qz
X = qt(X = q⊤zq⊤t)⊤· · · 形状最適化· · · 板厚分布最適化· · · 同時最適化(9)

fini = 16.506[kNm]
V = 202.98[m3 ]
図1解析モデル

3.2. 板厚分布・位相最適化
X = qt ，q̄t,min = −3m，q̄t,max = 5m，t¯max = 0.5m，α =500，(10)Vmin =0.8V̄ ，V = 202.98m3 として問題 (9) を解いて得られた板厚分布最適解を図 2, 図 3, 図 4 に示す．図中の f－ 520 －として表した値は，最適解におけるひずみエネルギーを示す．

3.3. 同時最適化結果()⊤X = q⊤z q⊤t ，q̄z,min = −30m，q̄z,max = 30m，q̄t,min =−3m，q̄t,max = 5m，t¯max = 0.5m，α =500，Vmin =0.8V̄ ，V = 202.98m3 として問題 (9) を解いて得られた最適化結

f = 1.376[kNm]
V = 162.381[m3 ]f / fini = 0.083
図2

果を図 6，図 7，図 8 に示す．

最適板厚分布 (CG 法)鳥瞰図f = 1.380[kNm]
V = 162.381[m3 ]f / fini = 0.084

立面図
図6

図3同時最適化結果 (CG 法)

最適板厚分布 (NewtonCG 法)

f = 1.376[kNm]
V = 162.381[m3 ]f / fini = 0.083
図4鳥瞰図

最適板厚分布 (加速勾配法)

立面図

Strain Energy[kNm]f = 1.210[kNm]
V = 162.371[m3 ]f / fini = 0.073

CGNewtonCG101agm

図7f = 0.887[kNm]
V = 162.380[m3 ]f / fini = 0.054

同時最適化結果 (NewtonCG 法)1000
図5255075100
Number of steps

125150175鳥瞰図

勾配法の違いによるひずみエネルギーと step 数の関係

板厚分布最適化結果 (図 2, 図 3, 図 4) から，支持点近傍の板厚を厚くし，支持点から最も離れた場所であるシェ

立面図
図8f = 0.890[kNm]
V = 162.381[m3 ]f / fini = 0.054

同時最適化結果 (加速勾配法)

ル中央部分の板厚が消失することで外力による曲げモーメントを小さくし，ひずみエネルギーを減少させていることが確認できる．ひずみエネルギーは初期形状と比べて 90% 以上減少した．勾配法の違いにかかわらず，概ね同様の板厚分布が得られた．また，ひずみエネルギーも概ね一致している．

報告 H95－ 521 －Strain Energy[kNm]CG
NewtonCG

101agm

今後の展望として，積載荷重の考慮が挙げられる．板厚が消失した範囲は自重は 0 となるが，積載荷重は一定値を取るためそのままでは変位が増大し構造的に不利となることから，積載荷重を考慮した状態で本手法をそのまま適用すると開口部を作り出すことができなくなると考えられる．今後は積載荷重も考慮可能な形に本手法を改

1000
図950100150
200
Number of steps

250300350

勾配法の違いによるひずみエネルギーと step 数の関係

同時最適化結果 (図 6，図 7，図 8) から，シェルのライズがカテナリー状に上がり，支持点近傍の板厚を厚く，支持点から最も離れた場所であるシェル中央部分の板厚が消失することでひずみエネルギーを減少させていることが確認できる．NewtonCG 法，加速勾配法を用いて全ポテンシャルエネルギー最小化問題を解くと，最適化結果において概ね同様のシェル形状，板厚分布，ひずみエネルギーが得られる．一方，CG 法を用いた場合，シェルのライズが低くひずみエネルギーが大きな解が得られた．ひずみエネルギーと step 数の関係 (図 5，図 9) から，全ポテンシャルエネルギー最小化問題に用いる勾配法によらず，収束性能に大きな差異が無いことが分かる．

4. 結言
本研究では，形状と板厚分布を B スプラインで離散化した上でシグモイド関数を板厚に乗じることで，連続最適化問題として定式化し，シェルの形状・板厚・位相最適化問題を解いた．また，剛性方程式を解く代わりに全ポテンシャルエネルギー最小化問題として定式化し，CG 法，NewtonCG 法，加速勾配法の 3 つの勾配法を適用し，目的関数性能と収束性能の比較を行った．研究結

良していきたい．
[参考文献]
1） 木村俊明, 大森博司. 形状と厚さの同時最適化法の定式化とその応用–自由曲面シェル構造の構造形態創生手法の提案 (その 1)–. 日
本建築学会構造系論文集, Vol. 74, No. 640, pp. 1091–1098, 2009.6.
2） 藤田慎之輔, 寒野善博. シグモイド関数を用いたシェル構造物の形状・板厚・位相同時最適化. 日本機械学会論文集, Vol. 85, No. 872,
pp. 18–00160, 2019.4.
3） Y.Kanno. Accelerated proximal gradient method for elastplastic analysis of large-scale trusses. Proceedings of OPTIS 2016 (Japan Societyof Mechanical Engineers), Vol. No.1213, , 2016.1.(in Japanese).
4） Y.Kanno S.Fujita and M.Ohsaki. Accelerated gradient method for geometrically nonlinear analysis of framed structures. Proceedings ofOPTIS 2016 (Japan Society of Mechanical Engineers), Vol. No.1212,, 2016.1.(in Japanese).
5） Y.Kanno. Accelerated proximal gradient method for equilibrium analysis of elastoplastic spatial truss structures. IASS Annual Symposium̶ Spatial Structures in the 21st Century, 2019.6.
6） R. E. Perez, P. W. Jansen, and J. R. R. A. Martins. pyopt: a pythonbased object-oriented framework for nonlinear constrained optimization. Structural and Multidisciplinary Optimization, Vol. 45, No. 1,pp. 101–118, 2012.1.
7） S.Fujita and Y.Kanno. Application of accelerated gradient method toequilibrium analysis of trusses with nonlinear elastic materials.j. struct.
constr. eng. AIJ, Vol. 84, No. 763.
8） Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland,Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson,Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew
Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew
R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhanPolat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, JosefPerktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R.
Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa,Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientiﬁc Computing in Python. Nature Methods,
Vol. 17, pp. 261–272, 2020.

果は次のようにまとめられる．

• 板厚にシグモイド関数を乗じたものを要素板厚とすることで，閾値以下の板厚は削除するという離散最適化問題を連続最適化問題に近似し，逐次二次計画法によって解くことができる．

• 全ポテンシャルエネルギー最小化問題として定式化し，節点変位ベクトルを求めることで，行列の特異化を防ぎ，概ね優良解が得られる．

• 板厚分布最適化の場合，勾配法の違いにかかわらず，概ね同様の板厚分布が得られる．

• 同時最適化の場合，全ポテンシャルエネルギー最小化問題に CG 法を用いると十分な収束性能がなく優良解が得られない．

報告 H95－ 522 －