法線方向の接触剛性

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

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

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

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

(64)

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

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

式（66）のマイナス符号は内力が外力のマイナスであることに由来します（外力とはマスター面からスレーブ節点に働く力です）。式（66）の法線をそのノルムで割ったヤコビベクトルで置き換え、u_i で微分（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)

これを微分すると以下の様になります（ξ と η は u_i の関数であること、 ${\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} }$ は未知で、後で求めます。これらは任意の u_i が変化したときの ξ と η の変化量を表します。k はスレーブ節点、またはマスター面に属する任意の節点です。ξ と η の値はスレーブ節点をマスター面に正射影して得られることに注意してください。

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

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

(73)

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

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

(74)

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

$\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}$ などが求まります。dq_i, i=1,..,n_m = 0 、dp_y = dp_z = 0 と仮定しましょう。すると上式の右辺は -q_{ξ_x}dp_x、-q_{η_x}dp_x となり、最終的に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 \frac{\partial \boldsymbol{F_j} }{\partial \boldsymbol{u_i} } = -... ...ight] - \varphi_j \frac{\partial \boldsymbol{F_p} }{\partial \boldsymbol{u_i} }$

(84)