FreeCADで片持ち梁の計算（静的応力解析）
今回はFreeCADでモデル化した【先端荷重を受ける片持ち梁】の静的応力解析例を紹介します
片持ち梁の形状寸法は断面幅B=50 mm，断面深さD=100 mm，梁長さL=1,000 mm です
材質はSS400相当とします
まずは梁理論を使って固定端の最大応力と先端のたわみを計算します
beam theory7.wxm
"Bernoulli-Eulerの梁理論"を用いた先端のたわみδEの算式を%o1に示します（式の導出については"梁理論 その4"を参照ください）
"Timoshenkoの梁理論"を用いた先端のたわみδTの算式を%o2に示します（式の導出については"梁理論 その5"を参照ください）
形状寸法 B=50 mm，D=100 mm，L=1,000 mm，ヤング率E=210,000 N/mm2，先端荷重P=10,000 Nをそれぞれ代入します（画面出力は省略）
断面2次モーメントIを計算した結果を%o8に示します（単位はmm4）
断面係数Zを計算した結果を%o9に示します（単位はmm3）
固定端の曲げモーメントMを計算した結果を%o10に示します（単位はN.mm）
最大垂直応力度σを計算した結果を%o11に示します（単位はN/mm2）
δEを計算した結果を%o12に示します（単位はmm）
ポアソン比ν = 0.3，せん断補正係数κ = 0.85をそれぞれ代入します（画面出力は省略）
せん断補正係数については"梁理論 その6"を参照ください
断面積Aを計算した結果を%o15に示します（単位はmm2）
横弾性係数Gを計算した結果を%o16に示します（単位はN/mm2）
δTを計算した結果を%o17に示します（単位はmm）
FreeCADのFEMワークベンチによる静的応力解析結果
メッシングにはGmsh，ソルバーにはCalculiXを使用しています
2Dのサーフェースモデル（beam_2d.FCStd）
Mises応力プロット（緑から赤にかけて応力が大きく，変位は10倍に強調表示しています）
Mises応力の最大値 119.47 N/mm2
Y方向変位の最大値 -3.82 mm
3Dのソリッドモデル（beam_3d.FCStd）
Mises応力プロット（緑から赤にかけて応力が大きく，変位は10倍に強調表示しています）
Mises応力の最大値 122.48 N/mm2
Y方向変位の最大値 -3.82 mm
ということで，FreeCADによる【先端荷重を受ける片持ち梁】の静的応力解析例を紹介しました
FreeCADはLGPLで作成されているフリーソフトなんですが，
モデリングからソルバー・プリポストまでワンパッケージで使えてしまうということに隔世の感がありますね(*´ω｀*)
ABC予想
今回はABC予想（abc conjecture）のお話です
2012年8月30日，京都大学の望月新一教授がABC予想を証明したとする論文を公開されました
この論文はABC予想の解決に宇宙際タイヒミュラー理論（IUT理論）が用いられており，今現在も査読中です
氏のblog”新一の「心の一票」”を読むとサーベイ作業は肯定的に終わっており，雑誌の最終査読報告書の提出が待たれます
IUT理論は革新的で大変難解な理論ですが，ABC予想自体はとても分り易いためmaximaで数値検証してみます
【abc予想】
a + b = c を満たす互いに素な自然数の組（a, b, c）（abc-triple）に対し，積 abc の根基(radical)を d と表す
このとき任意の ε > 0 に対して c > d^(1+ε) を満たす組（a, b, c）（例外的abc-triple）は高々有限個しか存在しない
abc conjecture.wxm
ABC予想の命題である不等式を%o1に示します
ここで根基は引数（自然数）の互いに異なる素因数の積を返す関数radとして%o2で定義します
ε= 0 として c < 10^4 の範囲で例外的abc-tripleがいくつあるか調べてみます
N：c の探索上限
abc-triple 15195742個 に対して例外的abc-tripleは 120個 となります
比率としては10^-5のorderであり十分小さいのですが，例外的abc-tripleを見てみると
a = 1 に対して，b = 8( = 3^2-1) や b = 80( = 3^4-1)，b = 6560( = 3^8-1) が見て取れます
ε= 0 の場合では (a = 1, b = 3^2^n-1, c = 3^2^n) のすべてのnについて例外的abc-tripleが無限に存在することが判っています
なのでABC予想は ε > 0 に対して述べられた命題となっています
ε= 0.4 として計算し直すと例外的abc-tripleは 120 → 3個まで減少します
ちなみに「ε= 1 とすると例外的abc-tripleは存在しない」という予想（強いabc予想）もあります
追記
2020年4月3日、京都大数理解析研究所の専門誌「PRIMS」に掲載されることが決まったとのことです
査読が通過するまで実に８年もの時間がかかってしまったことからも証明の難解さが伺えますね
巡回群
以前群論について群 その1や群 その2で簡単に触れましたが，今回は巡回群（cyclic group）について触れてみたいとおもいます
方程式 x^n = 1 の解を複素平面上にプロットすると正n角形が作図されます
（これは正17角形の作図可能性で扱ったのと全く同じ内容です）
cyclic group.wxm
いま n に5を代入します（%o1）（画面出力は省略）
%o2式に示す5次方程式の解を考えます
上式をxについて解いた結果を%o3式に示します
これらの解（1の複素5乗根）を複素平面上にプロットした結果を%t7に示します
解の一つを適当に選んでgとします（%o8）
このgを元の方程式の解集合Gに順次掛けていきます（G0→G1→G2→G3→G4）
ここでgは複素平面における回転変換に相当し，G5 = G0 となり巡回することが判ります
解集合Gの1番目の解（g = %e^(2*%i*%pi/5)）について考えます（%t11）
複素平面上の1に対してgを1回掛けた結果は原点まわりに2*%pi/5[rad]回転した赤点となります
これに順次gを掛けていくと青線を辿って巡回していき，5回目で元の1に戻ります
5乗すると1になるのはgが方程式x^5 = 1の解なので当然ですね
余談ですが，この"5"をgの位数(order)と呼びます
同様にして，
2番目の解（g = %e^(4*%i*%pi/5)）は4*%pi/5[rad]の回転（%t12）
3番目の解（g = %e^-(4*%i*%pi/5)）は-4*%pi/5[rad]の回転（%t13）
4番目の解（g = %e^-(2*%i*%pi/5)）は-2*%pi/5[rad]の回転（%t14）
5番目の解（g = 1）は0(= 2*%pi)[rad]の回転（%t15）を，それぞれ表します
ということで，1の複素5乗根が乗法に関して巡回群を成すことを示しました
追記
複素平面におけるこのような巡回群はU(1)と呼ばれます
通常の実平面における同様の巡回群はSO(2)と呼ばれます
SO(2)とU(1)は同型（isomorphism）です
吊荷重
久し振りの更新です（・ω・）
今回は吊荷重（Lifting Load）について簡単に説明します
lifting load.wxm
N : ワイヤー本数
M : 吊荷重量(kg)
ベクトル同士の外積と，ベクトルのnormを返す関数を%o4, %o5式でそれぞれ定義します（画面出力は省略）
ワイヤー8本1点吊りの問題を考えます(%o6)
吊荷として 4m x 2m x 100mmの鉄板(?)の重量をMに代入します(%o7)
吊荷の重心を原点としてワイヤー取付箇所8点分の位置ベクトル(m)を%o8〜15で代入します
H : 吊り点位置ベクトル(m)(%o16)
L : ワイヤー長さの配列(m)(%o17)
ワイヤー荷重ベクトルPを，荷重ベクトルの大きさを表すpを用いて%o18で定義します（単位はどちらもN）
%o19でパッケージ"draw"をロードします（画面出力は省略）
吊り点とワイヤーの吊り姿をdraw3dでプロットして確認します(%o20)
真上から見下ろした図を示しますが，実際には3D表示なのでグルグルと回せます
（複数点吊りする場合はN, H, P等を適当に修正してください）
W : 吊荷の荷重ベクトル(%o21)
pp : pの配列(%o22)
この問題の並進方向の釣合い方程式（equilibrium equation）をee1に代入します(%o23)
同様にして回転方向の釣合い方程式をee2に代入します(%o24)
ee1とee2を合わせてppについて解いた結果を%o25に示します
これより未知数の数に対して独立した方程式の数が足りず，助変数%r1〜%r3を3つ含んだ解が返されています
ワイヤー本数が多いほどこの助変数は増えていきます(´･ω･｀)
X : 助変数の配列(%o26)
%rnum_listは直前の結果内の助変数を配列として返してくれるとても便利な関数です
鉛直方向の荷重の2乗の残差をRとして定義します(%o27)
RについてXのヤコビヤンJを%o28で計算します（画面出力は省略）
Rの停留条件（J = 0）をXについて解いた結果を%o29に示します
助変数がすべて求まったので，ppを再評価します(%o30)
pの総和p[sum]はW[3]よりも大きな値となることに注意してください
一般に吊り角度が大きくなる（吊り高さH[3]が低くなる）とワイヤー荷重は大きくなります
得られた解でee1およびee2を再評価した結果，どちらの釣合い方程式も満足していることが解ります(%o32,33)
節点自由度の縮約 その2
前回の節点自由度の縮約を使って”半剛接接合部を有する梁要素”の導出について説明したいと思います
上図のように梁の両端に半剛接回転バネ（@で表示）が二つ接続した構造を考えます
- バネ要素①の曲げ剛性をKa，要素長さを0とします
- バネ要素③の曲げ剛性をKb，要素長さを0とします
- 梁要素②のヤング率をＥ，断面2次モーメントをI，要素長さをLとします
- 梁要素内で材料・形状ともに一定とします
節点1および4の自由度は鉛直方向変位vと回転角θの2つ
内節点2および3の鉛直方向は拘束されているので回転角θのみとなり，全自由度は 2 x 2 + 2 x 1 = 6 となります
condensation2.wxm
節点変位ベクトルを返す関数を%o1式で定義します
バネ要素節点変位ベクトルを返す関数を%o2式で定義します
要素①〜③の要素剛性マトリックスを%o4〜6式にそれぞれ示します
要素①〜③の要素節点力ベクトルを%o7, 9, 10式にそれぞれ示します
全体節点変位ベクトルUを%o11式に示します
節点力を全て重ね合わせてUで係数を括りだしたものが全体剛性マトリックスKとなります(%o12)
これが縮約していない通常の（全自由度 6 x 6 の）全体剛性マトリックスです
Uc : 縮約する節点変位ベクトル(%o14)
節点2と3の自由度の縮約を行います
まずはK.Uを計算します(%o13)
節点2と3に対応する自由度は3番目と4番目の成分です
節点2と3には外力が作用しないとして釣合い方程式を立て，Ucについて解いた結果が%o15式となります
%o15式の解を用いてK.Uを計算し直したものが%o16式です
3番目と4番目の成分がちゃんと 0 になっていることがわかります
また1，2，5，6番目の成分から縮約した変位成分θ2およびθ3が消えていることがわかります
Uから縮約される変位成分を除いた，節点1&4の節点変位ベクトルU14を%o17式に示します
残った節点力をU14で係数を括りだしたものが縮約された半剛接梁の要素剛性マトリックスKsとなります(%o18)
（全自由度は 6 - 2 = 4 となり，4 x 4 のマトリックスとなります）
式の簡単化のためにバネ要素の曲げ剛性を半剛接パラメータλ( = 0〜1 )を用いた%o21, 22式でそれぞれ定義します
半剛接パラメータλa, λbを用いてKsを書き直した結果を%o23式に示します
Ksに対してλa = 1, λb = 1 と仮定すると，両端剛接とした通常の要素剛性マトリックスと一致します(%o24)
Ksに対してλa = 1/2, λb = 1/2とすると，中間的な回転剛性を持つ要素剛性マトリックスを表現します(%o25)
Ksに対してλa = 0, λb = 0とすると，両端ピン接合としたトラス要素の剛性マトリックスと一致します(%o26)
追記
縮約された節点2および3の回転量を知りたい場合は別途計算する必要があります
%o15式の解をU14で係数を括りだしたものをNマトリックスとして%o27式に示します
半剛接パラメータλa，λbを用いて書き換えたものが%o28式となります
これより相対回転角θ' = θ14 - θ23 = θ14 - N.U14 で計算出来ます
ここで，θ14 = [θ1,θ4]，θ23 = [θ2,θ3]
節点自由度の縮約 その1
今回は節点自由度の縮約（Condensation）について説明したいと思います
上図のように梁の曲げ要素が二つ接続した構造を考えます
- 要素①の断面2次モーメントをI1、要素長さをr*Lとします
- 要素②の断面2次モーメントをI2、全体長さをLとします
- ヤング率はどちらの要素もＥで共通とします
- 要素内で材料・形状ともに一定とします
1節点の自由度は鉛直方向変位vと回転角θの2つ
節点数は合わせて3つなので全自由度は 2 x 3 = 6 となります
condensation1.wxm
曲げ梁の要素剛性マトリックスを返す関数kを%o1式で定義します
要素節点変位ベクトルを返す関数uを%o2式で定義します
要素①及び②の節点力ベクトルを計算した結果を%o3, 4にそれぞれ示します
全体節点変位ベクトルUを%o5式に示します
節点力を全て重ね合わせてUで係数を括りだしたものが全体剛性マトリックスKとなります(%o6)
これが縮約していない通常の（全自由度6x6の）全体剛性マトリックスです
Uc : 縮約する節点変位ベクトル(%o8)
節点2の自由度の縮約を行います
まずはK.Uを計算します(%o7)
節点2に対応する自由度は3番目と4番目の成分です
節点2には外力が作用しないとして釣合い方程式を立て，Ucについて解いた結果が%o9式となります
%o9式の解を用いてK.Uを計算し直したものが%o10式です
3番目と4番目の成分がちゃんと 0 になっていることがわかります
また1，2，5，6番目の成分から縮約した変位成分v2およびθ2が消えていることがわかります
Uから縮約される変位成分を除いた，節点1&3の節点変位ベクトルU13を%o11式に示します
式の簡単化のために無次元量Rを%o13式で定義します
残った節点力をU13で係数を括りだしたものが縮約された全体剛性マトリックスKcとなります(%o15)
（全自由度は 6 - 2 = 4 となり，4 x 4 のマトリックスとなります）
Kcに対して r = 1 と仮定すると，要素①の要素剛性マトリックスと一致します(%o16)
Kcに対して r = 0 と仮定すると，要素②の要素剛性マトリックスと一致します(%o17)
Kcに対して r = 1/2 と仮定すると，要素①と②が同じ長さ(L/2)で接続された場合の中間的な曲げ剛性を表します(%o18)
ルジャンドル多項式
今回はルジャンドル多項式（Legendre polynomial）のお話です
ルジャンドルの微分方程式の一般解に対して”λが非負整数かつ[1, 1]に確定点を持つ”という条件を課すことで与えられます
Gauss積分の所でこの多項式を内挿関数に使っていますね
legendre.wxm
n次のルジャンドル多項式P(n)は二項係数を用いて%o1式で表されます
%i1にてパッケージ"interpol"をロードします（画面出力は省略）
legendre_p関数を使って P(1)〜P(5) を計算した結果を示します(%t4〜8)
%t4〜8式を x = -1〜1の範囲でプロットします(%t9)
これより [1, 1]に確定点を持つことが解ります→P(n)|(x=1) = 1
ルジャンドル多項式が持つ性質の一つとして，”閉区間[−1, 1]上のL2-内積に関して直交する”というものがあります
%o10式に示す内積の積分を 10 x 10 まで具体的に計算した結果が%o11式となります
これより%o10式が直交性を持つこと，対角成分は正規化されておらず 2/(2*n+1) となることが判ります
