有限要素法の多点拘束の要素段階での組み込み

質問

FEMプログラムの中に多点拘束の機能を組み込みたいのですが、全体マトリクスの段階で組み込むと、従属関係を表すマトリクスが大きくなるので、要素マトリクスの段階で組み込む方法があれば、教えてくれないでしょうか。なお私の理解している全体マトリクスの場合は次です。
(1)元の剛性方程式：[K]{d}={F}
(2)従属関係：{d}=[C]{d0}、ここに[C]は従属関係を表すマトリクス、
{d0}は主（master)変位
ここで(2)を(1)に代入し、両辺に前から[C]の転値を乗じると
(3)　[C]T[K][C]{d0}=[C]T{F}
（３）式を解いて主変位を求め、さらに（２）から従変位を求める。

回答No.2

h191224 · ベストアンサー

回答が遅くなりました。

その「ある方」とは、この文献の著者ですね？
http://www.jsce.or.jp/committee/struct/civil_edu_home/prog/anketo/H19_femdoc02.pdf

MPCの従属変数を消去する関係式は、あなたの(3)式を見てもわかるように、あたかも座標変換の処理と同じようなマトリクスになります。
座標変換は、要素剛性マトリクスの作成の最中に処理するのが普通で、全体剛性マトリクスを作ってから変換処理をするのではありません。
ここでの処理は、はるかに複雑なので、「ある方」がなぜ「要素レベルで、等価な操作が可能」とおっしゃるのかわかりません。
その「ある方」は知り合いですので、今度お会いした時に、真意を確かめてみたいと思いますが、あなたへの回答には、とても間に合いませんね。

私が記載した方法(仮想要素法とでも呼びましょう)は、1970年代半ばに発表されましたが、 現在出版されている文献には、どこにも書かれていないと思います。
伊達政宗が、もう10年早く生まれていたら、天下をとったかも知れない、と言われるほどの逸材であったように、この仮想要素法も、素晴らしい方法なのですが、出てくるのが遅すぎました。
ですから、これを採用している老舗のプログラムは存在しないと思います。

私はこの仮想要素法が好きで、昔、自分でプログラムを作成していた頃は、好んで使用していました。 メモリーが少なかった頃は、この方法によると、メモリーを有効に使うことがで きて、処理が非常に楽だったのです。

ところで、あなたの追加された例では、節点４～６がｘｙの２自由度とも拘束されているので、全体剛性マトリクスは節点１～３に関するものだけが残り、次のようになります。（等幅フォントでご覧下さい。）

K11xx k11xy K12xx K12xy K13xx K13xy
K11yx k11yy K12yx K12yy K13yx K13yy
K21xx k21xy K22xx K22xy K23xx K23xy
K21yx k21yy K22yx K22yy K23yx K23yy
K31xx k31xy K32xx K32xy K33xx K33xy
K31yx k31yy K32yx K32yy K33yx K33yy

ここに
u1x-u2x=0
u2x-u3x=0
に対応する仮想要素をはめ込むと、マトリクスは次のようになります。

K11xx k11xy K12xx K12xy 1 K13xx K13xy 0
K11yx k11yy K12yx K12yy 0 K13yx K13yy 0
K21xx k21xy K22xx K22xy -1 K23xx K23xy 1
K21yx k21yy K22yx K22yy 0 K23yx K23yy 0
1 0 -1 0 0 0 0 0
K31xx k31xy K32xx K32xy 0 K33xx K33xy 1
K31yx k31yy K32yx K32yy 0 K33yx K33yy 0
0 0 1 0 0 -1 0 0

要は、
u1x-u2x=0に関係する自由度1x,2xが出現し終わった時点で、自由度1x,2xに対応する仮想要素をはめ込む、
u2x-u3x=0に関係する自由度1x,2xが出現し終わった時点で、自由度1x,2xに対応する仮想要素をはめ込む、
という操作を行えば良いのです。

この位置（以降）であれば、対角項が０であっても、何ら問題は生じません。
ですから、次のように、最後尾にまとめてはめ込んでも良いのです。

K11xx k11xy K12xx K12xy K13xx K13xy 1 0
K11yx k11yy K12yx K12yy K13yx K13yy 0 0
K21xx k21xy K22xx K22xy K23xx K23xy -1 1
K21yx k21yy K22yx K22yy K23yx K23yy 0 0
K31xx k31xy K32xx K32xy K33xx K33xy 0 -1
K31yx k31yy K32yx K32yy K33yx K33yy 0 0
1 0 -1 0 0 0 0 0
0 0 1 0 -1 0 0 0

お礼

再度のご丁寧な回答、ありがとうございます。
・「ある方の文献」、まさしくおっしゃる通りです。著書もある方なので、そんなうまい方法があるのかと思ってしまいました。
・ご回答の仮想要素法、ざっと見た感じではすばらしい方法のように思いますが、これから手で追いかけて正しい理解に努めたいと思います。

回答No.1

h191224

結論から言って、「要素マトリクスの段階で組み込む方法」というのは、実用上は無理です。

MPC(MultiPointConstraint＝多点拘束)を、要素段階で、ある要素のみに対して、組み込めるためには、関連する節点自由度が、ある要素内にのみ所属するということになります。
すると、その節点は、他の要素に所属できませんから、次の２つのいずれかに限定されます。
・他の領域と結合しない境界上の節点
・要素境界上にはなく、内部にのみ所属する節点

両方とも、極めて特殊なケースなのですが、本当にそれでよろしいのでしょうか？

Yesなら、話は簡単です。通常の従属節点自由度消去法と一緒です。
この場合には、全体剛性マトリクスの自由度は、その従属節点自由度分だけ減ります。

しかし、一般にMPCというものは、複数の要素に所属する節点自由度間で定義されます。
したがって、関連節点が複数要素にまたがるために、ある要素の段階でMPCを組み込んでしまうということは不可能です。要素剛性マトリクスを順番に作って行き、MPCに含まれる節点がすべて登場し終えるまで、MPCの方程式は処理し切れません。結果として、「従属関係を表すマトリクスが大きくなる」のと同じ現象が発生してしまいます。

なお、ご自分でプログラムを作成されていらっしゃるようなので、MPCの処理方法を１つご紹介しておきましょう。
従属・独立を設定せずに、全体剛性方程式に組み込んでしまう方法です。

具体例として、
ａｘ＋ｂｙ＋ｃｚ＝０
という３自由度間のＭＰＣを与えたいとしましょう。
この時に、仮の要素を作ってしまうのです。
その仮の要素の剛性マトリクスは、４×４の成分を持ち、その成分は次のように配置します。

０，０，０，ａ
０，０，０，ｂ
０，０，０，ｃ
ａ，ｂ，ｃ，０

余計な自由度が入ってきますが、従属自由度を消去して組み込む場合に対して、処理は楽で、全体剛性マトリクスもあまり増えません。
バンド幅は、従属自由度消去法と同様に増えますが、縦型消去法(注)を使っている場合には、全体剛性マトリクスの大きさや、処理速度には、ほとんど影響を与えず、さらに節点処理順にも影響を与えない点で優れています。
(3)では、[Ｃ]Tを乗じたためにマトリクスの非０成分が増えてしまいますが、この方法ではそのようなことも起きません。非常に良い方法です。

ただし、要注意点もあります。
ひとつは、ピボットが０なので、必ず上のように、ｎ行ｎ列目に、ピボットが来るようにして組み込む必要があります。
もうひとつは、同じMPCを２個以上与えても、気づかない点です。当然エラーになるので、データ作成の際の自己管理が必要です。

(注)通常の消去法は、行に沿って処理が進んでいきます。これが横型消去法です。
この横型消去法を少し変えるだけで、列に沿って処理を進めることができるようになるのです。これが縦型消去法です。
上記のようなタイプのマトリクスを全体剛性マトリクスに組み込むと、バンド幅が大きくなります。
しかし、このマトリクスは、ある成分から上がすべて０という性質を持っていて、組み込んだ後もその性質は変わりません。
縦型消去法は、このようなマトリクスの処理に非常に強みを発揮します。

なお、以上で、マトリクスの”要素”を、わざと”成分”と言い換えています。FEMの場合、混乱が起きるからです。

お礼

補足で、***以降の記述は勘違いでしたので、図以外は無視してください。
〔下も間違いです〕
[節点1－2，1－3]に対して
1,-1
-1,1

補足

ご丁寧な回答、ありがとうございます。私も「要素マトリクスの段階で組み込む方というのは、実用上は無理」とは思ったのですが、ネットである方が「要素レベルで、等価な操作が可能」と書いてあったので、もしやそういう方法があるのかと思ったしだいです。ご教示の方法、完全に理解したとはいいがたいのですが、単純な例では次ようになるのでしょうか（図がわかりにくくてすみません）。
・平面応力要素4個、自由節点：１～6、点2にx方向の力
・節点1,2,3のx方向の変位が等しいと仮定。
・仮要素剛性1(節点1，2)と仮要素剛性2(節点1，3)は同じ
・3行あるいは3列の成分は追加自由度に関する部分
・要素剛性
0,0,1
0,0,-1
1,-1,0
・上記要素剛性を全体剛性に重ね合わせる。
・なお、ピボットに関する記述がよくわかりません。また、この方法では対角がゼロになるため、連立方程式の求解でピボット選択をしなければいけないのでしょうか。
・この方法に関する参考文献を教えていただければありがたいのですが。
********************************************************
・これを書いていてふと思ったのですが、この例の場合、追加の自由度なしの次のような仮要素剛性ではだめでしょうか。これでもx1=x2=x3を満たすような気がするのですが。
[節点1－2，1－3]に対して
1,-1
-1,1

ｙ
↑　　⇒　 　　　　　
1・・・・・２・・・・・3
|------|------|
|------|------| 6
|//////|//////|→ｘ