法線方向の接触剛性
次へ 上へ 前へ 目次へ
次へ:Tangent contact stiffness 上へ:Node-to-Face Penalty Contact 前へ:General considerations 目次へ


法線方向の接触剛性

節点・面接触要素はマスター面につながったスレーブ節点で構成されます(図110を参照)。従ってこの接触要素は 1+nm 個の節点で構成されます。ここで nm はマスター面に属する節点の数です。有限要素の剛性行列は各節点での変位については各節点の内力の導関数になります。従って、私たちは節点での内力を決定する必要があります。

スレーブ節点位置を p、そのマスター面への投影位置を q とすると両者を結んだベクトルは以下のようなります。

r = p - q (61)

またこの位置での間隔 r は次のようになります。

r = rn (62)

ここで n はマスター面上のローカルな法線を意味します。マスター面に属する節点を$ \boldsymbol{q_i},i=1,n_m$、面内でのローカル座標を ξ 、ηとすると、次のように記述できます。

$\displaystyle \boldsymbol{q}= \sum _j \varphi _j(\xi ,\eta) \boldsymbol{q_j},$ (63)

$\displaystyle \boldsymbol{m} = \frac{\partial \boldsymbol{q} }{\partial \xi } \times \frac{\partial \boldsymbol{q} }{\partial \eta }$ (64)

さらに次の様になります。

$\displaystyle \boldsymbol{n} = \frac{\boldsymbol{m} }{\Vert \boldsymbol{m} \Vert}.$ (65)

m は面のヤコビベクトルです。節点 p での内力は以下の式に従います。

$\displaystyle \boldsymbol{F_p}= -f(r) \boldsymbol{n} a_p,$ (66)

ここで f はユーザーによって選択された圧力と間隔の関数、ap は節点 p に代表されるスレーブ領域面積です。N に属するスレーブ節点が面積 Ai のスレーブ面 i に触れる場合、その面積は以下の様に求められます。

$\displaystyle a_p=\sum _{i=1} ^N A_i/{n_s}_i.$ (67)

式(66)のマイナス符号は内力が外力のマイナスであることに由来します(外力とはマスター面からスレーブ節点に働く力です)。式(66)の法線をそのノルムで割ったヤコビベクトルで置き換え、ui で微分(i はスレーブ節点またはマスター面に属する任意節点)すると以下が得られます。

$\displaystyle \frac{1}{a_p} \frac{\partial \boldsymbol{F_p} }{\partial \boldsym...
...\otimes \frac{\partial \Vert \boldsymbol{m} \Vert}{\partial \boldsymbol{u_i} }.$ (68)

ここで

$\displaystyle \frac{\partial }{\partial \boldsymbol{u_i} } \left (\frac{\bolds...
...l{m}\Vert } \cdot \frac{\partial \boldsymbol{r} }{\partial \boldsymbol{u_i} } ,$ (69)

なので、上記の式は以下の様に書き換えられます。

$\displaystyle \frac{1}{a_p} \frac{\partial \boldsymbol{F_p} }{\partial \boldsymbol{u_i} } =$ $\displaystyle - \left(\frac{\partial f}{\partial r } \frac{1}{\Vert \boldsymbo...
...\frac{\partial \Vert \boldsymbol{m} \Vert }{\partial \boldsymbol{u_i} } \right]$
$\displaystyle + \frac{f}{\Vert \boldsymbol{m} \Vert } \left[ \boldsymbol{n} \ot...
...{u_i} } - \frac{\partial \boldsymbol{m} }{\partial \boldsymbol{u_i} } \right] .$ (70)

従って未定の導関数は$ \partial \boldsymbol{m}/\partial \boldsymbol{u_i}$$ \partial \boldsymbol{r}/ \partial \boldsymbol{u_i}$$ \partial \Vert \boldsymbol{m} \Vert / \partial \boldsymbol{u_i}$ということになります。

m の導関数は式(64)から得られ、以下の様に書くことができます。

$\displaystyle \boldsymbol{m}= \sum_j \sum_k \frac{\partial \varphi _j}{\partial...
...partial \varphi_k}{\partial \eta } [ \boldsymbol{q_j} \times \boldsymbol{q_k}].$ (71)

これを微分すると以下の様になります(ξ と η は ui の関数であること、$ {\partial \boldsymbol{q_i} }/{\partial \boldsymbol{u_j} }= \delta_{ij}\boldsymbol{I}$であることに注意してください)。

$\displaystyle \frac{\partial \boldsymbol{m} }{\partial \boldsymbol{u_i} } =$ $\displaystyle \left [ \frac{\partial^2 \boldsymbol{q} }{\partial \xi^2} \times ...
...partial \eta } \right] \otimes \frac{\partial \xi }{\partial \boldsymbol{u_i} }$
+ $\displaystyle \left[ \frac{\partial \boldsymbol{q} }{\partial \xi } \times \fra...
...partial \eta} \right] \otimes \frac{\partial \eta }{\partial \boldsymbol{u_i} }$
+ $\displaystyle \sum_{j=1}^{n_m} \sum_{k=1}^{n_m} \left[ \frac{\partial \varphi_j...
...right] (\boldsymbol{I} \times \boldsymbol{q_k}) \delta_{ij}, \:\:\: i=1,..n_m;p$ (72)

右辺にある導関数$ {\partial \xi }/{\partial \boldsymbol{u_i} }$$ {\partial \eta }/{\partial \boldsymbol{u_i} }$は未知で、後で求めます。これらは任意の ui が変化したときの ξ と η の変化量を表します。k はスレーブ節点、またはマスター面に属する任意の節点です。ξ と η の値はスレーブ節点をマスター面に正射影して得られることに注意してください。

式(61)と式(63)を組み合わせて r を得ると、ui の導関数は以下の様に書けます。

$\displaystyle \frac{\partial \boldsymbol{r} }{\partial \boldsymbol{u_i} }= \del...
...rtial \boldsymbol{u_i} } + \varphi_i (1 - \delta _{ip}) \boldsymbol{I} \right],$ (73)

ここで p はスレーブ節点を表します。

最終的にベクトルのノルムの導関数はベクトル自身の導関数で書き表すことができます。

$\displaystyle \frac{\partial \Vert \boldsymbol{m} \Vert }{\partial \boldsymbol{...
...ol{m} \Vert} \cdot \frac{\partial \boldsymbol{m} }{\partial \boldsymbol{u_i} }.$ (74)

残るは ξ と η の ui に対する導関数です。q がマスター面に対する p の正射影であるということは接続ベクトル r がマスター面に接するベクトル$ \partial \boldsymbol{q}/ \partial \xi $$ \partial \boldsymbol{q}/ \partial \eta $と直行するということです。

ここで

$\displaystyle \boldsymbol{r} \perp \frac{\boldsymbol{\partial q} }{\partial \xi }$ (75)

は以下の様に書くことができます。

$\displaystyle \boldsymbol{r} \cdot \frac{\partial \boldsymbol{q} }{\partial \xi } =0$ (76)

あるいは次の様に書くこともできます。

$\displaystyle \left [\boldsymbol{p} - \sum_i \varphi_i (\xi ,\eta) \boldsymbol...
...ft[ \sum_i \frac{\partial \varphi_i}{\partial \xi } \boldsymbol{q_i} \right]=0.$ (77)

上式を微分すると次の様になります。

$\displaystyle \left[ d \boldsymbol{p} - \sum_i \left(\frac{\partial \varphi_i}...
...ta + \varphi_i d \boldsymbol{q_i} \right) \right] \cdot \boldsymbol{q_{\xi }} +$
$\displaystyle \boldsymbol{r} \cdot \left[ \sum_i \left(\frac{\partial^2 \varph...
...\frac{\partial \varphi_i}{\partial \xi } d \boldsymbol{ q_i} \right) \right] =0$ (78)

ここで qξq の ξ に対する微分です。上式は以下の様に書き直すことができます。

$\displaystyle (d \boldsymbol{p} - \boldsymbol{q}_{\xi } d \xi - \boldsymbol{q}...
... } d \eta - \sum_i \varphi_i d \boldsymbol{q_i}) \cdot \boldsymbol{q}_{\xi } +$
$\displaystyle \boldsymbol{r} \cdot (\boldsymbol{q}_{\xi \xi} d \xi + \boldsymb...
...\eta + \sum_i \frac{\partial \varphi_i}{\partial \xi } d \boldsymbol{q_i})=0 .$ (79)

最終的に以下の式が得られます。

$\displaystyle (- \boldsymbol{q}_{\xi} \cdot \boldsymbol{q}_{\xi } + \boldsymbo...
...boldsymbol{q}_{\xi} + \boldsymbol{r} \cdot \boldsymbol{q}_{\xi \eta }) d \eta =$
$\displaystyle - \boldsymbol{q}_{\xi} \cdot d \boldsymbol{p} + \sum_i \left[ (\...
...tial \varphi_i}{\partial \xi } \boldsymbol{r}) \cdot d \boldsymbol{q_i} \right]$ (80)

η-方向の接線についても同様です。

$\displaystyle (- \boldsymbol{q}_{\xi} \cdot \boldsymbol{q}_{\eta } + \boldsymb...
...ldsymbol{q}_{\eta} + \boldsymbol{r} \cdot \boldsymbol{q}_{\eta \eta }) d \eta =$
$\displaystyle - \boldsymbol{q}_{\eta} \cdot d \boldsymbol{p} + \sum_i \left[ (...
...ial \varphi_i}{\partial \eta } \boldsymbol{r}) \cdot d \boldsymbol{q_i} \right]$ (81)

これから、$ \partial \xi / \partial \boldsymbol{q_i}$$ \partial \xi / \partial \boldsymbol{p}$などが求まります。dqi, i=1,..,nm = 0 、dpy = dpz = 0 と仮定しましょう。すると上式の右辺は -qξxdpx、-qηxdpx となり、最終的に2つの未知数$ \partial \xi / \partial p_x$$ \partial \eta / \partial p_x$を持つ2つの方程式になります。$ \partial \xi / \partial \boldsymbol{p}$が決まれば以下の式から自動的に$ \partial \xi / \partial \boldsymbol{u_p}$が決定します。

$\displaystyle \frac{\partial \xi }{\partial \boldsymbol{p} } = \frac{\partial \xi }{\partial \boldsymbol{u_p} },$ (82)

他の導関数についても同様です。これによって導関数$ \partial \boldsymbol{F}_p / \partial \boldsymbol{u_i}$が決定されます。

以下の式

$\displaystyle \boldsymbol{F_j}= - \varphi_j (\xi, \eta) \boldsymbol{F_p},$ (83)

から、次が得られます。

$\displaystyle \frac{\partial \boldsymbol{F_j} }{\partial \boldsymbol{u_i} } = -...
...ight] - \varphi_j \frac{\partial \boldsymbol{F_p} }{\partial \boldsymbol{u_i} }$ (84)

これがマスター節点での力の導関数になります。


次へ 上へ 前へ 目次へ
次へ:Tangent contact stiffness 上へ:Node-to-Face Penalty Contact 前へ:General considerations 目次へ
guido dhondt 2016-03-08