打撃試験における空洞形状判別に対する
レベルセット関数を用いたトポロジー最適化解析の適用
Application of topology optimization analysis using level set function for judgement of cavity shape in hammering test
正 倉橋 貴彦
*1，村上 祐貴
*2，正 外山 茂浩
*2， 正 池田 富士雄
*2，正 井山 徹郎
*2，正 井原 郁夫
*1 Takahiko KURAHASHI*1, Yuki MURAKAMI*2, Shigehiro TOYAMA*2,Fujio IKEDA*2, Tetsuro IYAMA*2 and Ikuo IHARA*1
*1 長岡技術科学大学 Nagaoka University of Technology
*2 長岡工業高等専門学校 National Institute of Technology, Nagaoka College
代表著者：倉橋貴彦
[email protected]キーワード：逆問題とデータ同化の最新展開, トポロジー最適化, レベルセット関数, 構造内の空洞, 非破壊 打撃試験
1. 序 論
2012年の山梨県笹子トンネル崩落事故を契機に，コンクリートの打音試験に対する定量的評価に関する研究が 近年も注目されている．2020年6月12日の日刊工業新聞の一面には，「高速道路端 管理に新手法 データサイ エンスの活用」という記事もあり，分析ツールを使用した研修会も実施されている．現在存在しているインフラ の構造物は 50 年を経過したものも増加傾向にあり，少子高齢化の問題も踏まえると計算力学の技術を活かした 非破壊検査ツールの構築が急務であると考えられる．
コンクリート内に剥離が生じた際，剥離の形状に依っては，応力集中によるき裂端部における応力値も異なる ことから，剥離の形状はコンクリート片の落下を左右することにもなる．そのため，空洞の形状を把握すること はコンクリートの剥離防止の観点からも大変重要なものとなる．打音検査時にコンクリート表面において得られ た応答波形を用いて，空洞形状を推定する逆問題を考えた場合，逆解析手法としては，形状最適化やトポロジー 最適化の技術を適用することができる．空洞の形状推定に形状最適化の解析法を用いる場合，逆解析開始時にメ ッシュの境界により空洞をあらかじめ表現しておき，その境界形状を反復計算により更新をすることにより，空 洞形状の推定を実施することになる(1)．そのため，あらかじめ指定した空洞形状の個数を増やすことができない ことから，形状の個数を逆解析前に何等かの方法で把握できておく必要がある．一方、トポロジー最適化では，
逆解析前に，形状を細かく表現できる分解能を持ったメッシュを事前に準備しておき，各要素における材料の有 無を反復計算により判別できるため，空洞の個数がいくつある場合においても，空洞の形状推定が可能となる．
トポロジー最適化の解析手法では，密度法やレベルセット関数を用いた方法(2)が比較的良く用いられている．密 度法はヤング率を密度パラメータで表現することにより，密度パラメータの分布を求める解析を行うことになる が，密度パラメータを更新するために用いる感度の値を直接用いた場合，空洞と材料を有する領域がチェッカー ボード状に得られるため，結果として多孔質な形状が得られることが知られている．この点を打破するために，
感度に対するフィルタリングを適用し密度パラメータの分布の更新が行われるが，感度分布が平滑化されるため，
空洞部分と材料を有する境界を分け完全に二値化することができず，空洞と材料部分がミックスされた曖昧な境 界形状を推定することになる．この手法に対しレベルセット関数を用いた方法では，レベルセット関数に関する 反応拡散方程式を解くことによりレベルセット関数の正負により空洞部の境界形状を表現するため，密度法に比 べて空洞の境界形状を明確に表現できる方法として知られている．
以上のことを踏まえ，本研究では，レベルセット関数を用いたトポロジー最適化に着目し，コンクリート内の 空洞の形状を推定する解析について数値実験結果を示す．また，表面における打撃により得られる変位の時間履 歴を測定値として用いる．測定値については，空洞形状を仮定したメッシュに対して振動解析をすることにより，
あらかじめ人工的に測定値を算定しておく．そして，異なる空洞形状から逆解析を実施し，最初に仮定した形状 の推定の可否について検討を行う．
2. レベルセット関数を用いたトポロジー最適化
まず，評価関数を定義する．本研究では，表面の応答変位が評価対象となるため，シミュレーションの結果と 測定値の残差二乗値を観測点において積分したものを評価関数とする．式(1)に示すように評価関数を定義し，u は変位ベクトル，uobs.は変位の測定値によるベクトルを示す．また，Qは重み対角行列であり，測定点では対角成 分に1，そうでなければ，対角成分に0を入力する．t0は初期時刻，tfは終端時刻を示す．
𝐽 =1
2∫ (𝒖 − 𝒖𝑡𝑡𝑓 𝒐𝒃𝒔.)𝑇𝑸(𝒖 − 𝒖𝒐𝒃𝒔.)𝑑𝑡
0 （1）
振動方程式により，変位ベクトルuが算定されるため，評価関数の制約条件として有限要素法により離散化した 振動方程式(3)を考慮すると式(2)に示すラグランジュ関数が得られる．ここに，ドットは時間微分を示し，λは随 伴変数ベクトルを示す．また，本研究では，レベルセット関数の正負により領域の有無を判定するため，振動方 程式を有限要素法により離散化した際，要素領域の積分により得られる係数行列M，C，Kをレベルセット関数 の関数とし，𝑴(𝜙)，𝑪(𝜙)，𝑲(𝜙)と表すことにする．fは外力ベクトルを示す．また，レベルセット関数の更新式 の説明も考え，式(2)の非積分関数をFとしておく．
𝐽∗= 𝐽 + ∫ 𝝀𝑡𝑡𝑓 𝑇(𝑴(𝜙)𝒖̈ + 𝑪(𝜙)𝒖̇ + 𝑲(𝜙)𝒖 − 𝒇)𝒅𝒕
0 = ∫ 𝐹𝑑𝑡𝑡𝑡𝑓
0 （2）
式(2)に示すラグランジュ関数J*に対する第一変分を計算し，停留条件δJ*=0を計算する．ラグランジュ関数の随 伴変数に関する微分より式(3)に示す振動方程式が得られ，またラグランジュ関数の各変位成分に対する微分より 式(4)に示す随伴方程式が得られる．式(3)，式(4)の時間方向の離散化にはニューマークのβ法を適用する．
𝑴(𝜙)𝒖̈ + 𝑪(𝜙)𝒖̇ + 𝑲(𝜙)𝒖 = 𝒇 （3）
𝑴𝑻(𝜙)𝝀̈ − 𝑪𝑻(𝜙)𝝀̇ + 𝑲𝑻(𝜙)𝝀 = −𝑸𝑻(𝒖 − 𝒖𝒐𝒃𝒔.) （4）
ここで，式(5)に示す反応拡散方程式によりレベルセット関数が進展するものと仮定する．式中においてカンマ は微分を表す．式(5)の離散化に有限要素法を適用し，要素領域により離散化すると，式(6)に示す有限要素方程式 が得られる．ここに，𝜅(𝜙)は正のパラメータ，Cは𝐹𝑒,𝜙の大きさを調節するためのパラメータ，τは正則化パラメ ータ，𝑡̃は最適化のための仮想的な時間，Neは形状関数ベクトルを示す．仮想時間𝑡̃による微分に対しては，後退 差分法を適用し離散化を行う．
𝜙, 𝑡̃− 𝜅(𝜙)𝜏(𝜙,𝑥𝑥+ 𝜙,𝑦𝑦+ 𝜙,𝑧𝑧) = −𝜅(𝜙)𝐶𝐹,𝜙 （5）
𝑴𝑒𝝓𝑒, 𝑡̃+ 𝜅(𝜙)𝜏𝑺𝑒𝝓𝑒= −𝜅(𝜙)𝐶𝐹𝑒,𝜙∫ 𝑵Ω 𝑒𝑑Ω
e （6）
式(5)および式(6)における𝐹𝑒,𝜙は上述のように，係数行列M，C，Kをレベルセット関数𝜙により微分することによ り計算する（式(7)）．ここでは，通常の係数行列M，C，Kにレベルセット関数𝜙が乗じられているものとし，式 (8)のように表す．ここに， 𝐻𝑒(𝜙)は各要素におけるヘビサイド関数を示す．
𝐹𝑒,𝜙 = 𝝀𝑒𝑇(𝑴𝑒,𝜙𝒖𝑒̈ + 𝑪𝑒,𝜙𝒖𝑒̇ + 𝑲𝑒,𝜙𝒖𝑒) （7）
𝐹𝑒,𝜙= 𝝀𝑒𝑇(𝑴𝑒𝒖𝑒̈ + 𝑪𝑒𝒖𝑒̇ + 𝑲𝑒𝒖𝑒)𝐻𝑒(𝜙) （8）
ここに計算の流れを以下に示す．
1. 計算メッシュ，境界条件，初期条件，計算パラメータを入力する．また反復回数をk=0と設定する．
2. 振動方程式（式(3)）を解いて，評価関数（式(1)）を計算する．
3. k>0のとき，評価関数の収束判定を行う．|J(k)－J(k-1)/J0|<εの場合，計算を終了し，そうでない場合は次のステ ップへ進む．
4. 随伴方程式（式(4)）を解いて，式(8)を各要素において計算する．
5. レベルセット関数の時間進展式（式(6)を全要素において重ね合わせた式）により，各節点におけるレベルセ ット関数を計算する．
6. レベルセット関数が正の場合は物体領域，負の場合は空洞領域，0ならば物体の境界とし，構造形態を更新す る．k=k+1とし，ステップ2に戻る．
3. 数値解析例
数値解析例として用いるモデル図を図1に示す．打撃点は上部中央の点とし，変位の測定点は4点設定した．
解析に用いた数値パラメータを表1に整理する．本検討では，図2の右の目標形状に示すように，空洞が打撃点 の左右に1つずつ，計2つの空洞が存在するものとし，トポロジー最適化解析における形状の初期設定は図2の 左のように設定した．打撃条件としては，ガウシアンパルスF(t)=Fmaxexp(-(t-tpeak)2/s2) により与えることとし，Fmax
は2000N，tpeakは10-3s，sは10-4sとする．変位の経時変化を図2の右の目標形状モデルに対して振動解析を行う
ことで得られる振動変位を人工的な測定値として用い，図2の左に示す初期形状から空洞形状の推定を行う．本 検討では正則化パラメータτをCase-A : τ=0.001，Case-B : τ=0.005，Case-C : τ=0.010と変え，推定された空洞形状 に関して比較を行う．振動変位の測定値は，各測定点においてz方向変位のみ適用する．
図3に評価関数の履歴を示す．τの値が小さい程，評価関数の値は高いところで収束しており，逆にτの値が大 きい程，評価関数の値を下げられる結果が得られている．各々のτの条件において，推定された空洞の形状を図 4に示す．τの値が小さい程，空洞の領域は広がる傾向にあり，τの値が大きい程，空洞の領域が小さくなる結果 が得られた．図2の右に示す目標形状と比べると，推定された結果においてCase-C : τ=0.010の時が空洞の位置と しては最も近い結果が得られたが，空洞の形状は目標形状に比べて小さく得られていることを確認できる．本検 討の結果より，τ の値は空洞の大きさを調節するパラメータとして設定はできるが，トポロジー最適化解析にお いて設定する数値パラメータや打撃点や変位応答の測定点の位置によっては，正確に空洞のサイズを推定するの が困難になるものと考えられる．
4. 結 論
本研究では，レベルセット関数を用いたトポロジー最適化解析を非破壊打撃試験に適用し，空洞形状の推定可 能性に関する検討を実施した．結果として，トポロジー最適化解析において用いる正則化パラメータτを調節す ることにより，空洞の形状が変わることを確認した．正則化パラメータτの値によっては，正解とする空洞位置 に近い結果を得ることもできたが，空洞のサイズについては，良好な一致を示す結果は得られなかった．本解析 による空洞形状の推定は，打撃試験により得られる表面応答の測定点の位置にも依存すると考えられるため，測 定点の位置を変えた場合の検討も必要となると考えられる．また，本研究では，測定値として，数値解析による
結果を人工的な測定値として用いているため，実際のコンクリート試験体に対する打撃試験により得られた変位 応答波形を測定値として使用し，実問題への適用可能性についても検討を実施する必要がある．
Fig.1 Setup of design domain and size of computational domain
Table 1 Numerical conditions
Number of nodes 112761
Number of elements 100000
Number of time steps 256
Time increments Δt 39.0625
Young’s modulus E, GPa 35.096
Poisson’s ratio ν 0.16
Mass density ρ, kg/m3 2300
Damping coefficient cM, cK ( Reyleigh damping ) 90.0, 1.0×10-6 Weight parameters qu, qv, qw ( Diagonal component in matrix Q ) 0, 0, 1.0×109
Virtual time increment Δ𝑡̃ 1.0×10-8
Convergence criterion ε 1.0×10-3
Initial shape Target shape Fig.2 Initial and target shapes of cavity
×
x y
z
Fig.3 Variation of normalized performance function
Case-A : τ=0.001 Case-B : τ=0.005 Case-C : τ=0.010 Fig.4 Comparison of cavity shape for each regularization parameter τ
謝 辞
本研究は新潟県建設技術センターならびに長岡技術科学大学学長戦略的経費による研究助成を受け実施した．
研究遂行に際して長岡技術科学大学大学院機械創造工学専攻2019年度修了生の吉原健太氏の多大な協力を得た．
数値解析においては九州大学情報基盤研究開発センターの高性能演算サーバーシステムを利用した．記して謝意 を表す．
文 献
(1) 吉原健太，倉橋貴彦，村上祐貴，外山茂浩，池田富士雄，井山徹郎，井原郁夫，“ランダム・トンネリング・アル ゴリズムを用いた打撃応答波形に基づく三次元空洞位置・形状同定解析”, 日本機械学会北陸信越支部 第57期総 会・講演会講演論文集，(2020), pp.1-5.
(2) 山田崇恭, 西脇眞二, 泉井一浩, 吉村允孝, 竹澤晃弘, “レベルセット法による形状表現を用いたフェーズフィール ド法の考え方に基づくトポロジー最適化”, 日本機械学会論文集A 編, Vol.75, (2009), pp.550-558.
(3) 竹内則雄, 樫山和男, 寺田賢二郎, 計算力学 －有限要素法の基礎－，森北出版株式会社, (2003).
0.00E+00 2.00E-01 4.00E-01 6.00E-01 8.00E-01 1.00E+00 1.20E+00
0 5 10 15 20 25
Value of normalized performance function
Number of iterations
Case-A : τ=0.001 Case-B : τ=0.005 Case-C : τ=0.010
X Y
Z
X Y
Z
X Y
Z
