Contact Problem with Kuhn-Tucher Condition

Kuhn-Tucher条件による簡単な接触問題

Last Update:2009年05月15日

概要

弾塑性体や形状最適化問題や接触問題など，制約条件が不等式で表される場合の最適化問題はKuhn-Tucker条件(KKT条件)を使うことができる．Kuhn-Tucker条件はラグランジュ未定乗数法の不等式への拡張である．

Kuhn-Tucher条件

次のような問題を考える．

$\min\;W(\b{u})\;while\;G(\b{u})\le 0$

このような解はKuhn-Tucker条件から以下を満たす

$\{\begin{array}{l}\frac{\partial W}{\partial \b{u}}+\mu \frac{\partial G}{\partial \b{u}}=0\\ \mu \ge 0\\ \mu G(\b{u}) = 0\end{array}$

ここでの $\mu$ が拘束力を表す． $\mu\ge 0$ は拘束力が常に負でないこと． $\mu G(\b{u}) = 0$ はか $\mu$ のどちらかが０であることを意味している．例えばある点が壁でさえぎられている場合は、拘束力は壁と反発する方向に壁の上にあるときにしか働かないことに対応している．

未定乗数の意味

関数の選び方は幾つも考えられるが

$|\frac{\partial G}{\partial \b{u}}| = 1$

を満たすようにしておくと，

$\frac{\partial W}{\partial \b{u}}+\mu\frac{\partial G}{\partial\b{u}}=0$

$|\frac{\partial W}{\partial\b{u}}|=|\mu||\frac{\partial G}{\partial \b{u}}|=\mu$

であるから， $\mu$ は制約条件の満たされない方向へのエネルギーの勾配の大きさに等しくなる．エネルギーの変位での微分は力であるから，この場合 $\mu$ は拘束力の大きさを意味している．このようにLagrange未定乗数に拘束力の大きさとしての意味づけをすることができるので便利である．

多数の不等式拘束条件

不等式拘束条件が沢山ある場合はそれぞれ別々に考えればよい．

$\min\;W(\b{u})\;while\;G^k(\b{u})\le 0\;\;\;k=\{1,\ldots,n\}$

このような解はKuhn-Tucker条件から以下を満たす

$\{\begin{array}{l}\frac{\partial W}{\partial \b{u}}+\mu^k \frac{\partial G^k}{\partial \b{u}}=0\\ \mu^k \ge 0\\ \mu^k G^k(\b{u}) = 0\end{array}$

有限要素法への組み込み

問題定式化

不等式拘束条件が課されている節点の数がだけあるとする．また，番目の拘束条件の節番号がであるとする．

$\min\;W(\b{u})\;while\;G(\b{u}^{a(k)})\le 0\;\;\;\{k=1,\ldots,n\}$

上の議論から有限要素法の解は以下を満たす

$\{\begin{array}{l}\frac{\partial W}{\partial \b{u}^a}+\mu^a \frac{\partial G}{\partial \b{u}^a}=0\\ \mu^a \ge 0\\ \mu^a G(\b{u}^a) = 0\end{array}$

この式が満たされるようにする解法を考えよう

反復解法

このような条件は非線形なので、反復計算を行い上の式が満たされるようにしなければならない．例えば次のような反復のさせ方が考えられる． $\delta,\epsilon$ を十分小さな数とする．

$\mu>\delta$ なら前回Lagrange未定乗数を入れて計算した結果，今回も $\mu>0$ である．だからLagrange未定乗数を入れる．
１以外で， $G>\epsilon$ なら， $G\le 0$ が満たされないので，になるように，Lagrange未定乗数を入れる．
１.２以外の時，つまり $\mu\le\delta$ かつ $G\le\epsilon$ の時は $G\le 0$ が満たされると考え，Lagrange未定乗数を入れずに、 $\mu$ を０にセットする

これを図にすると図１のようになる．青い場所は接触が起こっているとして、 G=0 が満たされるようなLagrange未定乗数を入れる場所である．赤い場所は変位について拘束を入れず $\mu$ に０をセットする場所である．最終的に黒い線の上に落ち着くことで、 $G\le 0,\;\mu\ge 0,\; G\mu =0$ が満たされるようになれば反復終了である．

図１．ＫＫＴ条件の反復解法の条件わけ図．青い部分は制約条件を考慮して，赤い部分は制約条件を考慮しない．

条件が不連続的に切り替わる部分は数式で書けば綺麗だが，プログラミング上は誤差とか考慮しなければならない．要は浮動小数点同士の比較は正確にできないので上の $\delta,\epsilon$ が必要である．この数字はパラーメータスタディによって決めなければならない．

このように2つの状況が不連続に存在する状況では、反復中に2つが交互に入れ替わって収束が遅くなるchatteringという現象が起こる可能性がある．

有限要素法定式化

上の反復法を有限要素法で定式化する時の例を示す．上の反復では $\mu$ との値によって２つの場合に分けられた．

次の増分に関する連立一次方程式を解き，解を更新する

$\[\begin{array}{ll}[K_{uu}]^{ab}_{ij}&[K_{u\mu}]^{ab}_{i}\\[K_{\mu u}]^{ab}_{\;j}&[K_{\mu\mu}]^{ab} \end{array}\] \{\begin{array}\{\Delta u\}^b_j\\\{\Delta\mu\}^b\end{array}\} = \{\begin{array}\{f\}^b_j\\0\end{array}\} - \{\begin{array}\{Q_u\}^b_j\\\{Q_\mu\}^b\end{array}\}$

$\{\begin{array}\b{u}^{new}=\b{u}^{old}+\Delta\b{u}\\\mu^{new}=\mu^{old}+\Delta\mu\end{array}$

但し， $\{Q_u\},\{Q_\mu\}$ をそれぞれ，変位と拘束力についての等価節点内力であるとした．

ラグランジュ未定乗数を入れる場合

$\mu>\delta$ もしくは $G>\epsilon$ の場合はラグランジュ未定乗数を入れて G=0 が満たされるようにする．

内力については次のようになる．

$\{\begin{array}{l}\{Q_u\}^a_i=\frac{\partial W}{\partial\b{u}^a_i}+\mu^a\frac{\partial G}{\partial\b{u}^a_i}\\\{Q_\mu\}^a=G(\b{u}^a)\end{array}$

これを微分することで接線剛性行列が得られる

$\{\begin{array}{l}[K_{uu}]^{ab}_{ij}=\frac{\partial \{Q_u\}^a_i}{\partial {u}^b_j} = \frac{\partial^2 W}{\partial\b{u}^a_i\partial\b{u}^b_j}+\mu^a\frac{\partial^2 G}{\partial\b{u}^a_i\partial\b{u}^b_j}\delta_{ab}\\ [K_{u\mu}]^{ab}_{i}=\frac{\partial \{Q_u\}^a_i}{\partial {\mu}^b} = \frac{\partial G}{\partial\b{u}^a_i}\delta_{ab}\\ [K_{\mu u}]^{ab}_{\;j}=\frac{\partial \{Q_\mu\}^a}{\partial u^b_j } = \frac{\partial G}{\partial\b{u}^b_j}\delta_{ab}\\ [K_{\mu \mu}]^{ab}=\frac{\partial \{Q_\mu\}^a}{\partial \mu^b } = 0 \end{array}$

ラグランジュ未定乗数を入れない場合

$\mu<\delta$ かつ $G<\epsilon$ の場合は，節点は拘束力を受けない．よって $\mu=0$ にする．

$\{\begin{array}{l}\{Q_u\}^a_i=\frac{\partial W}{\partial\b{u}^a_i}\\\{Q_\mu\}^a=\mu^a\end{array}$

これを微分することで接線剛性行列が得られる

$\{\begin{array}{l}[K_{uu}]^{ab}_{ij}=\frac{\partial \{Q_u\}^a_i}{\partial {u}^b_j} = \frac{\partial^2 W}{\partial\b{u}^a_i\partial\b{u}^b_j}\\ [K_{u\mu}]^{ab}_{i}=\frac{\partial \{Q_u\}^a_i}{\partial {\mu}^b} = 0\\ [K_{\mu u}]^{ab}_{\;j}=\frac{\partial \{Q_\mu\}^a}{\partial u^b_j } = 0\\ [K_{\mu \mu}]^{ab}=\frac{\partial \{Q_\mu\}^a}{\partial \mu^b } = \delta_{ab} \end{array}$

直線との接触解析

直線と弾性体との接触を考えてみた．「ある直線よりも節点が上にある」という制約条件は不等式で書き表すことができる．

例えば点 $\b{p}$ を通り，単位法線ベクトルが $\b{n}$ である直線より上にあるという条件は

$\;\;G(\b{u})=-(\b{x}-\b{p})\cdot\b{n}\le 0\\ \Leftright G(\b{u})= -(\b{X}+\b{u}-\b{p})\cdot\b{n}\le 0$

と書ける．但し， $\b{x}$ は変位後の点の座標で， $\b{X}$ は点の初期座標． $\b{u}$ は初期状態からの変位を表す．但し，の微分は次のようになる．

$\frac{\partial G(\b{u})}{\partial u_i} = -n_i$

$\frac{\partial^2 G(\b{u})}{\partial u_i\partial u_j} = 0$

線形弾性体と直線との接触解析例

実際にこれらを実装してみたのが以下の動画である．

剛体円との接触

剛体円との接触を考える.円の外側に節点があると考えればよい．は円の半径で， $\b{p}$ は円の中心の座標とすると，拘束条件の関数は次のようになる．

$\;\;\;G(\b{u})=d-|\b{x}-\b{p}|\le 0\\ \Leftright G(\b{u})=d-|\b{X}+\b{u}-\b{p}|\le 0$

と書ける．但し， $\b{x}$ は変位後の点の座標で， $\b{X}$ は点の初期座標． $\b{u}$ は初期状態からの変位を表す．

$\frac{\partial G(\b{u})}{\partial u_i} = -\frac{(\b{u}+\b{X}-\b{p})_i}{|\b{u}+\b{X}-\b{p}|}$

$\frac{\partial^2 G(\b{u})}{\partial u_i\partial u_j} = -\frac{\delta_{ij}}{|\b{u}+\b{X}-\b{p}|}+\frac{(\b{u}+\b{X}-\b{p})_i(\b{u}+\b{X}-\b{p})_j}{|\b{u}+\b{X}-\b{p}|^3}$

が非線形性を持っており，２回微分が０にならないことに注意されたい．

線形弾性体と剛体円の接触解析例

下向きの体積力をかけている線形弾性体の梁に対して左端にリサージュ曲線の強制変位を与えて円と接触させている．

St.Venant-Kirchhoff体と剛体円の接触解析例

St.Venant-Kirchhoff体でできた梁を自由落下させて剛体円と接触させてみる，有限要素法解析をDelFEMに組み込んで実験してみた．

梁が円に当たってボヨーンと脇に跳ね返った後，回転しながら落ちてゆく．

St.Venant-Kirchhoff体は幾何学的非線形性を考慮した構成式で書かれるために，線形弾性体と違って回転しても不自然な膨張が起こらない．

参考文献

Wikipedia Karush-Kuhn-Tucker_conditions: http://en.wikipedia.org/wiki/Karush-Kuhn-Tucker_conditions

クーン・タッカーの必要条件: http://www.isigas.com/KuhnTucker.html

最適化2006 Kuhn-Tucker条件: http://www.fbc.keio.ac.jp/~hkomiya/education/lecture/saitekika-2006-3.pdf

Kuhn-Tucker条件を利用した多変数関数における最適化: http://dspace.wul.waseda.ac.jp/dspace/bitstream/2065/13120/1/Honbun-t3606u116.pdf