摩擦ブリック要素との有限すべり/有限変形接触

変形可能な構造の有限すべり/有限変形接触に対処するための基本的なツールはGWFL（汎用弱形式言語）でアクセス可能です．いくつかの補間変換（ 補間変換 を参照）は接触検出を実行し，反対側の境界境界にある接触境界から積分できるように定義されています．実際の構成における単位法線ベクトルやCoulomb摩擦との接触を考慮する投影などの他の有用なツールも弱形式言語の演算子として定義されています．

もちろん，有限すべり/有限変形接触アルゴリズムの計算コストは微小すべり-微小変形よりも大幅に高くなります．

レイトレーシング補間変換

高水準汎用構築に接触検出を組み込むために，特定の補間変換が定義されています（内挿変換の詳細については 補間変換 を参照してください）． [KO-RE2014] に記載されているレイトレーシング接触検出に基づいており，以下で説明する基準を使用します．補間変換は，異なる可能性のある接触面を記憶します．ほとんどの方法で，潜在的接触面はマスターとスレーブの2つのカテゴリに分類されます（ 図 を参照）．

スレーブ表面は "コンタクタ" で，マスターは "ターゲット" です．剛支障も考慮されます．それらは常にマスター表面です．基本的な規則は，スレーブ表面とマスタ面の間で接触が考慮されることです．しかし，マルチコンタクトフレームオブジェクトと GetFEM 要素は2つのマスター表面間の接触，マスター表面と任意の数のスレーブ表面とマスター表面の自己接触を含む複数接触の状況を可能にします．

基本的には，接触対を検出するために，Gauss点またはスレーブ表面の有限要素法節点はマスター表面に投影されます（ 図 を参照）．自己接触が考慮される場合，Gauss点またはマスター表面の有限要素法節点もマスター表面に投影されます．

レイトレーシング変換のモデルへの追加は次の通りです．

void add_raytracing_transformation(model &md, const std::string &transname,
                                    scalar_type d)

ここで transname はGWFLでそれを参照できる変換に与えられた名前で， d は解放距離です（上記参照）．

レイトレーシング変換はスレーブまたはマスターの接触境界なしで追加されます．次の関数は変換にいくつかの境界を追加することを可能にします．

add_master_contact_boundary_to_raytracing_transformation(model &md,
           const std::string &transname, const mesh &m,
           const std::string &dispname, size_type region)

add_slave_contact_boundary_to_raytracing_transformation(model &md,
           const std::string &transname, const mesh &m,
           const std::string &dispname, size_type region)

where dispname is the variable name which represents the displacement on that contact boundary. The difference between master and slave contact boundary is that the contact detection is to be performed starting from a slave or master boundary toward a master boundary. The contact detection is not performed toward a slave boundary. Consequently, only the influence boxes of the elements of the master surfaces are computed and stored.

次の関数で剛支障（マスター表面とみなされる）を追加することも可能です．

add_rigid_obstacle_to_raytracing_transformation(model &md,
           const std::string &transname,
           const std::string &expr, size_type N)

where expr is the expression of a signed distance to the obstacle using the syntax of GWFL (X being the current position, X(0), X(1) ... the corresponding components). For instance an expression X(0) + 5 will correspond to a flat obstacle lying on the right of the position -5 of the first coordinate. Be aware that the expression has to be close to a signed distance, which in particular means that the gradient norm has to be close to 1.

In order to distinguish between non-contact situations and the occurrence of a contact with another deformable body or with a rigid obstacle, the transformation returns an integer identifier which can be used by the Interpolate_filter command of GWFL (see 補間変換). The different values:

• 0：このGauss点で接点が見つかりません

• 1：このGauss点で変形可能な物体と接触します

• 2：このGauss点で剛体の障害物と接触します

これらの3つの場合，次のように微分することが可能です．

Interpolate_filter(transname, expr1, 0)
Interpolate_filter(transname, expr2, 1)
Interpolate_filter(transname, expr3, 2)

GWFLでは， expr1 ， expr2 ， expr3 は計算されるべき異なる項に対応しています． matlabインタフェースのデモプログラム /interface/tests/matlab/demo_large_sliding_contact.m は，使用例を示しています．

modelオブジェクトが使用されていない場合，変換は ga_workspace オブジェクトで直接使用することもできます．詳細は， getfem/getfem_contact_and_friction_common.h を参照してください．また，modelオブジェクトのフレームワークでは，この変換のインターフェース使用が，以下に説明するモデルブリック要素によって許可されることにも注意してください．

接触対検出アルゴリズム

接触対は，スレーブ（または自己接触の場合はマスター）の点と，最も近いマスター表面（または剛性の障害物）の投影点とによって形成されます．使用されたアルゴリズムは 図 のように整理されます．

（[Pantz2008] など，)グローバルなトポロジカルな基準なしに有効な接触先と無効な接触先の間で確実に微分することは不可能です．しかし，この種の基準は，実施するのに非常にコストがかかる可能性があります．したがって，一般的には，可能性のあるすべてのケースをカバーすることができない簡単なヒューリスティックな基準を実装します．このような基準をここに提示します．これらはもちろん不完全であり，変更の対象となります．最初に， 図 では，条件を微分しなければならない有効なまたは無効な接触先の数の状況を見ることができます．

アルゴリズムの詳細は次の通りです．

• 影響ボックスの計算 要素の影響ボックスは，放出距離に等しい距離でその境界ボックスのオフセットだけです．この戦略を使用する場合，リリース距離は要素サイズに比べて大きすぎるべきではありません．さもなければ，点は接触対の検索をかなり減速させる多数の影響ボックスに対応します．インフルエンスボックスは，領域ツリーオブジェクトに格納され，平均複雑度を有するアルゴリズムを有する点を含むボックスを， \(O(log(N))\) で見つけます．

• 潜在的な接触対とは何か． 潜在的な接触対は調査されるスレーブ点 - マスターエレメント面のペアです．マスタ表面上のスレーブ点の投影が行われ，基準が適用されます．

• 投影アルゴリズム． マスター要素面へのスレーブ点の投影は，幾何変換と変位場による参照要素上の表面のパラメータ化によって行われます．投影中，要素面内にとどまる制約は適用されません．これは，要素面が解析的に延長されることを意味します．投影は，パラメータ化およびNewtonおよび/またはBFGSアルゴリズムを使用して，スレーブ点と投影点との間の距離を最小にすることによって実行されます． raytrace がtrueに設定されている場合，投影は計算されません．その代わりに，点xからの単位法線ベクトルの方向にyを見つけるためにレイトレーシングします．これは，通常の状況の逆を意味します（xはyの投影です）．

基準のリストは

• 基準1：ユニットの法線円錐/ベクトルは互換性があり，2つの点は同じ要素を共有しません． 2つの単位法線ベクトルはスカラー積が非正であれば互換性があります． f.e.mの場合．フェムト節点は一般に複数の要素で共有されているため，各要素の単位法線ベクトルからなる法線円が考慮されます．少なくとも1組の単位法線ベクトルがスカラ積が非正であれば，2つのノーマルコーンは互換性があります．計算を単純化するために，ノーマルコーンの立体角がマルチコンタクトフレームオブジェクトのパラメータ cut_angle よりも小さい場合，ノーマルコーンは平均法線ベクトルに縮小されます．この基準は，問題（B）および（K1）に対応することを可能にします．

• 基準2：投影/レイトレース点の検索が収束しない場合，接触対は削除されます． マスターエレメント上のスレーブ点の投影/レイトレースを計算するために使用されたNewton法（および投影用のBFGS）表面が収束しない場合，対は考慮されず警告が生成されます．

• 基準3：投影された点は要素の内側に設ける． スレーブ点は，面の内側に留まるという制約なしにマスター要素の表面に投影されます（つまり，面が長くなることを意味します）．正射影が顔の外側にある場合，ペアは考慮されません．しかし，現状では，問題（J3）に対応するためには，随意の対応が考慮されなければなりません（この点では拘束を受けて顔に投影し，この時点で正常な円錐を試験する）．この基準は， （F2），（K2），（M1），（M2）のケースに対応できます．

• 基準4：リリース距離が適用されます． スレーブ点とマスター表面上の投影間の距離がリリース距離より大きい場合，接触対は考慮されません．リリース距離が適応され，変形がそれほど重要でない場合，（C），（E），（F1），（G），（H）の場合を処理できます．

• 判定基準5：堅い障害物との比較 スレーブ点とマスタ表面上の投影との間の符号付き距離が，堅い障害物を有するものより大きい場合（リリース距離が剛性障害物にも最初に適用されることを考慮して）接触対は考慮されません．

• 基準6：自己接触のみ：基準配置における単位法線の試験を適用します． 自己接触の場合，スレーブ点とマスター要素が同じメッシュに属し，スレーブ点はマスタ面の背後にあり（ユニットの外向きの法線ベクトルに対して），リリース距離の4倍も離れていません．これは，問題（A），（C），（D），（H）に対応することができます．

• 基準7：接触対の最小符号付距離 保持された接触対（または剛性障害物）の間で，符号付き最小距離に対応するものが保持されます．

投影による節点接触要素

表記： \(\Omega \subset \rm I\hspace{-0.15em}R^d\) は，変形可能な物体の基準配置を示し，いくつかの未接続部分（ 図 を参照）で構成されています．\(\Omega_t\) は変形後の構成で \(\varphi^h: \Omega \rightarrow \Omega_t\) は有限要素法空間 \(V^h\) の変形近似です．変形 \(u^h: \Omega \rightarrow \rm I\hspace{-0.15em}R^d\) は \(\varphi^h(X) = X + u^h(X)\) により定義されます．基準配置 \(\Omega\) の一般的な点は \(X\) によって表記されており， \(x = \varphi^h(X)\) によって変形後の対応する点が表記されいます．変形後の構成の対応する境界はそれぞれ \(\Gamma_t^S\) と \(\Gamma_t^M\) です．(変形後構成の)境界での点 \(x = \varphi^h(X)\) における外向き単位基準化ベクトルは \(n_x\) により表現されています．最後に，表記 \(\delta A[B]\) は変形に関する量 \(A\) と方向 \(B\) の指向性導関数を表す．同様に，表記 \(\delta^2 A[B,C]\) は方向 \(B\) と方向 \(C\) の2階導関数を示す．

\(J(\varphi^h)\) は，接触と摩擦の寄与を考慮しない，システムのポテンシャルエネルギーです．典型的には，それは弾性及び外部負荷ポテンシャルエネルギーを含みます．\(X_i\) は次のようにします． \(i \in I_{\text{nodes}}\) は参照設定のスレーブ境界上の有限要素節点の集合です． \(i \in I_{\text{def}}\) は，変形可能な物体のマスター表面と接触する可能性のある接触節点です． \(X_i\) は次のようにします． \(i \in I_{\text{rig}}\) は，堅い障害物と接触する可能性のある接触節点です．

\(x_i = \varphi^h(X_i)\) は，変形された構成の対応する節点であり， \(y_i\) は変形された構成のマスター表面（または剛体の障害物）の投影です． \(Y_i\) はマスター表面上点を確認する \(y_i = \varphi^h(Y_i)\) です．これにより，通常のギャップを次のように定義できます．

\[g_i = n_y . (\varphi^h(X_i) - \varphi^h(Y_i)) = \|\varphi^h(X_i) - \varphi^h(Y_i)\| \text{Sign}(n_y . (\varphi^h(X_i) - \varphi^h(Y_i))),\]

ここで， \(n_y\) は， \(y\) のマスター表面の外向きの単位法線ベクトルです．

定常剛性障害物のみを考慮し，Alart-Curnier Lagrangian [AL-CU1991] の原理を適用すると，摩擦条件との節点接触の問題は非対称バージョンで次のように表すことができる（線形弾性のケース）．

\[\begin{split}\left\{\begin{array}{l} \mbox{Find } \varphi^h \in V^h \mbox{ such that } \\ \displaystyle \delta J(\varphi^h)[\delta u^h] - \sum_{i \in I_{\text{def}}} \lambda_i \cdot (\delta u^h(X_i) - \delta u^h(Y_i)) - \sum_{i \in I_{\text{rig}}} \lambda_i \delta u^h(X_i) = 0 ~~~ \forall \delta u^h \in V^h, \\ \displaystyle \dfrac{1}{r} \left[\lambda_i + P_{n_y, {\mathscr F}}(\lambda_i + r\left(g_i n_y - \alpha(\varphi^h(X_i) - \varphi^h(Y_i) - W_T(X_i)+W_T(Y_i)))\right)\right]= 0 ~~\forall i \in I_{\text{def}}, \\[1em] \displaystyle \dfrac{1}{r} \left[\lambda_i + P_{n_y, {\mathscr F}}(\lambda_i + r\left(g_i n_y - \alpha(\varphi^h(X_i) - W_T(X_i)))\right)\right]= 0 ~~\forall i \in I_{\text{rig}}, \end{array}\right.\end{split}\]

ここで， \(W_T, \alpha, P_{n_y, {\mathscr F}}\) ... + 接線系

申し訳ありませんが，このブリック要素は動作していません．

摩擦接触のための高水準汎用構築のツール

次の非線形演算子は，GWFLで定義されています（ 任意の項を計算する - 高水準の汎用的な構築手順 - 弱形式言語 (GWFL) 参照）．

• Transformed_unit_vector(Grad_u, n) ここで Grad_u は変位フィールドの勾配であり， n は基準配置の単位ベクトルです．この非線形演算子は，以下の通りです．

  \[n_{trans} = \dfrac{(I+ \nabla u)^{-T} n}{\|(I+\nabla u)^{-T} n\|}\]

  偏微分は以下の通りです．

  \[ \begin{align}\begin{aligned}\partial_{u} n_{trans}[\delta u] = -(I - n_{trans}\otimes n_{trans})(I+ \nabla u)^{-T}(\nabla \delta u)^T n_{trans}\\\partial_{n} n_{trans}[\delta n] = \dfrac{(I+ \nabla u)^{-T}\delta n - n_{trans}(n_{trans}\cdot \delta n)}{\|(I+\nabla u)^{-T} n\|}\end{aligned}\end{align} \]

• Coulomb_friction_coupled_projection(lambda, n, Vs, g, f, r) ここで lambda は接触力であり， n は単位法線ベクトルであり， Vs はすべり速度であり， g はギャップ， f は摩擦係数， r は正の増強パラメータです．演算子の式は次のとおりです．

  \[ \begin{align}\begin{aligned}P(\lambda, n, V_s, g, f, r) = -(\lambda\cdot n + rg)_- n + P_{B(n,\tau)}(\lambda - rV_s)\\\mbox{with } \tau = \mbox{min}(f_3 + f_1(\lambda\cdot n + rg)_-, f_2)\end{aligned}\end{align} \]

  ここで， \((\cdot)_-\) は負の部分 (\((x)_- = (-x)_+\)) です:math:f_1, f_2, f_3 は摩擦係数3成分です．成分 \(f_2, f_3\) はオプションであることに注意してください．仮想係数が（ \(f_1\) のみ）与えられている場合，これは古典Coulomb摩擦則に相当します． 2つの成分（ \(f_1, f_2\) のみ）ベクトルが与えられている場合，これは与えられた閾値とのCoulomb摩擦に対応します．最後に，3つの成分からなるベクトルが与えられた場合，摩擦法則は，上で与えられた \(\tau\) の式に対応します．

  The expression \(P_{B(n,\tau)}(q)\) refers to the orthogonal projection (this is linked to the return mapping algorithm) on the tangential ball (with respect to \(n\) of radius \(\tau\).

  微分は \(T_n = (I - n \otimes n)\) および \(q_{_T} = T_n q\) のように表すことができます．

  \[ \begin{align}\begin{aligned}\begin{split}\partial_q P_{B(n,\tau)}(q) = \left\{\begin{array}{cl} 0 & \mbox{for } \tau \le 0 \\ \mathbf{T}_n & \mbox{for } \|q_{_T}\| \le \tau \\ \dfrac{\tau}{\|q_{_T}\|} \left(\mathbf{T}_n - \dfrac{q_{_T}}{\|q_{_T}\|}\otimes \dfrac{q_{_T}}{\|q_{_T}\|} \right) & \mbox{otherwise } \end{array} \right.\end{split}\\\begin{split}\partial_{\tau} P_{B(n,\tau)}(q) = \left\{\begin{array}{cl} 0 & \mbox{for } \tau \le 0 \mbox{ or } \|q_{_T}\| \le \tau \\ \dfrac{q_{_T}}{\|q_{_T}\|} & \mbox{otherwise} \end{array} \right.\end{split}\\\begin{split}\partial_n P_{B(n,\tau)}(q) = \left\{ \begin{array}{cl} 0 & \mbox{for } \tau \le 0 \\ -q \cdot n~\mathbf{T}_n - n \otimes q_{_T} & \mbox{for } \|q_{_T}\| \le \tau \\ -\dfrac{\tau}{\|q_{_T}\|} \left( q \cdot n \left(\mathbf{T}_n - \dfrac{q_{_T}}{\|q_{_T}\|}\otimes \dfrac{q_{_T}}{\|q_{_T}\|} \right) + n \otimes q_{_T} \right) & \mbox{otherwise.} \end{array} \right.\end{split}\\\partial_{\lambda} P(\lambda, n, V_s, g, f, r) = \partial_q P_{B(n,\tau)} +\partial_{\tau}P_{B(n,\tau)} \otimes \partial_{\lambda} \tau +H(-\lambda\cdot n - r\,g)~n \otimes n,\\\begin{split}\partial_{n} P(\lambda, n, V_s, g, f, r) = \left|\begin{array}{l} \partial_n P_{B(n,\tau)} +\partial_{\tau} P_{B(n,\tau)} \otimes \partial_n \tau \\ \hspace*{3em}+H(-\lambda\cdot n - r\,g) ~ \left(n \otimes \lambda - (2~\lambda\cdot n + r\,g)~n \otimes n + (\lambda\cdot n + r\,g)~\mathbf{I}\right), \end{array}\right.\end{split}\\\partial_{g} P(\lambda, n, V_s, g, f, r) = \partial_{\tau} P_{B(n,\tau)} ~ \partial_g \tau +H(-\lambda\cdot n - r\,g)~r~n\\\partial_{f} P(\lambda, n, V_s, g, f, r) = \partial_{\tau} P_{B(n,\tau)} \partial_{f} \tau\\\partial_{r} P(\lambda, n, V_s, g, f, r) = H(-\lambda\cdot n - r\,g)gn + \partial_q P_{B(n,\tau)}V_s +\partial_{\tau} P_{B(n,\tau)} \partial_r \tau\end{aligned}\end{align} \]

レイトレース付き積分接触要素

ブリック要素を追加します．

indbrick = add_integral_large_sliding_contact_brick_raytracing
  (model &md, const std::string &dataname_r,
   scalar_type release_distance,
   const std::string &dataname_friction_coeff = "0",
   const std::string &dataname_alpha = "1");

このブリック要素は，複数の接触先の状況に対処できます．前のセクションで説明したように，モデルにレイトレーシング補間変換を追加します．その名前はコマンドによって取得できます．

const std::string &transformation_name_of_large_sliding_contact_brick(model &md,
                       size_type indbrick);

ブリック要素がモデルに追加されると，マスターとスレーブの接触境界を次の関数で追加する必要があります．

add_contact_boundary_to_large_sliding_contact_brick(model &md,
    size_type indbrick, const mesh_im &mim, size_type region,
    bool is_master, bool is_slave, const std::string &u,
    const std::string &lambda = "", const std::string &w = "",
    bool frame_indifferent = false)

ここで， region は境界を表す有効なメッシュ領域番号でなければなりません． is_master はその接触境界で接触検出が行われるなら true にセットされ， is_slave はその境界上で接触項の積分が行われる場合は true に設定してください．接触境界は，特に自己接触検出を可能にするために，マスターとスレーブの両方になることが可能であることに留意してください． u は変位変数です． is_slave がtrueに設定されている場合， lambda は接触境界上の自由度を持つ乗数変数を記述しなければなりません（通常は md.add_filtered_fem_variable(...) メソッドでモデルに追加されます）．純粋なマスター接触境界は，乗数の定義を必要としません．さらに， w は進展する場合であり，前の時間ステップにおける変位を表します．

剛性の障害物を要素に加えることができます:

add_rigid_obstacle_to_large_sliding_contact_brick(model &md,
    size_type indbrick, std::string expr, size_type N)

ここで， expr はGWFL（ X は現在の位置）を使った式で，障害物までの符号付き距離でなければなりません． N はメッシュ寸法です．