Hatena::ブログ(Diary)

Maximaでこうぞうりきがく

2018-08-26

吊荷重

久し振りの更新です(・ω・)

今回は吊荷重(Lifting Load)について簡単に説明します

lifting load.wxm

f:id:ryooji_f:20180826090334p:image

N : ワイヤー本数

M : 吊荷重量(kg)

ベクトル同士の外積と,ベクトルのnormを返す関数を%o4, %o5式でそれぞれ定義します(画面出力は省略)

ワイヤー8本1点吊りの問題を考えます(%o6)

吊荷として 4m x 2m x 100mmの鉄板(?)の重量をMに代入します(%o7)

吊荷の重心を原点としてワイヤー取付箇所8点分の位置ベクトル(m)を%o8〜15で代入します


https://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20180826/20180826090331.png

H : 吊り点位置ベクトル(m)(%o16)

L : ワイヤー長さの配列(m)(%o17)

ワイヤー荷重ベクトルPを,荷重ベクトルの大きさを表すpを用いて%o18で定義します(単位はどちらもN)


f:id:ryooji_f:20180826090329p:image

f:id:ryooji_f:20180826090326p:image:w500

%o19でパッケージ"draw"をロードします(画面出力は省略)

吊り点とワイヤーの吊り姿をdraw3dでプロットして確認します(%o20)

真上から見下ろした図を示しますが,実際には3D表示なのでグルグルと回せます

(複数点吊りする場合はN, H, P等を適当に修正してください)


https://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20180826/20180826090323.png

W : 吊荷の荷重ベクトル(%o21)

pp : pの配列(%o22)

この問題の並進方向の釣合い方程式(equilibrium equation)をee1に代入します(%o23)

同様にして回転方向の釣合い方程式をee2に代入します(%o24)

ee1とee2を合わせてppについて解いた結果を%o25に示します

これより未知数の数に対して独立した方程式の数が足りず,助変数%r1〜%r3を3つ含んだ解が返されています

ワイヤー本数が多いほどこの助変数は増えていきます(´・ω・`)


https://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20180826/20180826090321.png

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]が低くなる)とワイヤー荷重は大きくなります


f:id:ryooji_f:20180826090318p:image

得られた解でee1およびee2を再評価した結果,どちらの釣合い方程式も満足していることが解ります(%o32,33)

2017-07-02

節点自由度の縮約 その2

前回の節点自由度の縮約を使って”半剛接接合部を有する梁要素”の導出について説明したいと思います


    Element
 Node1[1]2      [2]       3[3]4
     *-@-*―――――――――――――*-@-*
         |__________________________|
                      L

上図のように梁の両端に半剛接回転バネ(@で表示)が二つ接続した構造を考えます

  • バネ要素,龍覆温篝をKa,要素長さを0とします
  • バネ要素の曲げ剛性をKb,要素長さを0とします
  • 梁要素△離筌鵐偉┐鬘邸っ婆2次モーメントをI,要素長さをLとします
  • 梁要素内で材料・形状ともに一定とします

節点1および4の自由度は鉛直方向変位vと回転角θの2つ

内節点2および3の鉛直方向は拘束されているので回転角θのみとなり,全自由度は 2 x 2 + 2 x 1 = 6 となります

condensation2.wxm

http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170702/20170702225338.png

節点変位ベクトルを返す関数を%o1式で定義します

バネ要素節点変位ベクトルを返す関数を%o2式で定義します

要素 銑の要素剛性マトリックスを%o4〜6式にそれぞれ示します


f:id:ryooji_f:20170702225331p:image

要素 銑の要素節点力ベクトルを%o7, 9, 10式にそれぞれ示します


f:id:ryooji_f:20170702225328p:image

全体節点変位ベクトルUを%o11式に示します


f:id:ryooji_f:20170702225320p:image

節点力を全て重ね合わせてUで係数を括りだしたものが全体剛性マトリックスKとなります(%o12)

これが縮約していない通常の(全自由度 6 x 6 の)全体剛性マトリックスです


http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170702/20170702225316.png

Uc : 縮約する節点変位ベクトル(%o14)

節点2と3の自由度の縮約を行います

まずはK.Uを計算します(%o13)

節点2と3に対応する自由度は3番目と4番目の成分です

節点2と3には外力が作用しないとして釣合い方程式を立て,Ucについて解いた結果が%o15式となります


http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170702/20170702225309.png

%o15式の解を用いてK.Uを計算し直したものが%o16式です

3番目と4番目の成分がちゃんと 0 になっていることがわかります

また1,2,5,6番目の成分から縮約した変位成分θ2およびθ3が消えていることがわかります

Uから縮約される変位成分を除いた,節点1&4の節点変位ベクトルU14を%o17式に示します


http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170702/20170702225306.png

残った節点力をU14で係数を括りだしたものが縮約された半剛接梁の要素剛性マトリックスKsとなります(%o18)

(全自由度は 6 - 2 = 4 となり,4 x 4 のマトリックスとなります)


http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170702/20170702225302.png

式の簡単化のためにバネ要素の曲げ剛性を半剛接パラメータλ( = 0〜1 )を用いた%o21, 22式でそれぞれ定義します

半剛接パラメータλa, λbを用いてKsを書き直した結果を%o23式に示します


f:id:ryooji_f:20170702225258p:image

Ksに対してλa = 1, λb = 1 と仮定すると,両端剛接とした通常の要素剛性マトリックスと一致します(%o24)

Ksに対してλa = 1/2, λb = 1/2とすると,中間的な回転剛性を持つ要素剛性マトリックスを表現します(%o25)

Ksに対してλa = 0, λb = 0とすると,両端ピン接合としたトラス要素の剛性マトリックスと一致します(%o26)


追記

http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170702/20170702225252.png

縮約された節点2および3の回転量を知りたい場合は別途計算する必要があります

%o15式の解をU14で係数を括りだしたものをNマトリックスとして%o27式に示します

半剛接パラメータλa,λbを用いて書き換えたものが%o28式となります

これより相対回転角θ' = θ14 - θ23 = θ14 - N.U14 で計算出来ます

ここで,θ14 = [θ1,θ4],θ23 = [θ2,θ3]

2017-06-27

節点自由度の縮約 その1

今回は節点自由度の縮約(Condensation)について説明したいと思います


            Element
   Node1      [1]      2     [2]      3
       *===============*―――――――*
       |_______________|       |
       |_____ r*L ____________________|
                       L

上図のように梁の曲げ要素が二つ接続した構造を考えます

  • 要素,涼婆2次モーメントをI1、要素長さをr*Lとします
  • 要素△涼婆2次モーメントをI2、全体長さをLとします
  • ヤング率はどちらの要素もEで共通とします
  • 要素内で材料・形状ともに一定とします

1節点の自由度は鉛直方向変位vと回転角θの2つ

節点数は合わせて3つなので全自由度は 2 x 3 = 6 となります

condensation1.wxm

http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170627/20170627225224.png

曲げ梁の要素剛性マトリックスを返す関数kを%o1式で定義します

要素節点変位ベクトルを返す関数uを%o2式で定義します

要素ゝ擇哭△寮疆昔ベクトルを計算した結果を%o3, 4にそれぞれ示します


f:id:ryooji_f:20170627225221p:image

全体節点変位ベクトルUを%o5式に示します


http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170627/20170627225217.png

節点力を全て重ね合わせてUで係数を括りだしたものが全体剛性マトリックスKとなります(%o6)

これが縮約していない通常の(全自由度6x6の)全体剛性マトリックスです


http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170627/20170627225212.png

Uc : 縮約する節点変位ベクトル(%o8)

節点2の自由度の縮約を行います

まずはK.Uを計算します(%o7)

節点2に対応する自由度は3番目と4番目の成分です

節点2には外力が作用しないとして釣合い方程式を立て,Ucについて解いた結果が%o9式となります


http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170627/20170627225202.png

%o9式の解を用いてK.Uを計算し直したものが%o10式です

3番目と4番目の成分がちゃんと 0 になっていることがわかります

また1,2,5,6番目の成分から縮約した変位成分v2およびθ2が消えていることがわかります

Uから縮約される変位成分を除いた,節点1&3の節点変位ベクトルU13を%o11式に示します


http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170627/20170627225157.png

式の簡単化のために無次元量Rを%o13式で定義します

残った節点力をU13で係数を括りだしたものが縮約された全体剛性マトリックスKcとなります(%o15)

(全自由度は 6 - 2 = 4 となり,4 x 4 のマトリックスとなります)


f:id:ryooji_f:20170627225152p:image

Kcに対して r = 1 と仮定すると,要素,陵彖嚢篝マトリックスと一致します(%o16)

Kcに対して r = 0 と仮定すると,要素△陵彖嚢篝マトリックスと一致します(%o17)

Kcに対して r = 1/2 と仮定すると,要素,鉢△同じ長さ(L/2)で接続された場合の中間的な曲げ剛性を表します(%o18)

2017-06-14

ルジャンドル多項式

今回はルジャンドル多項式(Legendre polynomial)のお話です

ルジャンドル微分方程式一般解に対して”λが非負整数かつ[1, 1]に確定点を持つ”という条件を課すことで与えられます

Gauss積分の所でこの多項式を内挿関数に使っていますね


legendre.wxm

f:id:ryooji_f:20170614214122p:image

n次のルジャンドル多項式P(n)は二項係数を用いて%o1式で表されます


f:id:ryooji_f:20170614214119p:image

%i1にてパッケージ"interpol"をロードします(画面出力は省略)

legendre_p関数を使って P(1)〜P(5) を計算した結果を示します(%t4〜8)


f:id:ryooji_f:20170614214114p:image

%t4〜8式を x = -1〜1の範囲でプロットします(%t9)

これより [1, 1]に確定点を持つことが解ります→P(n)|(x=1) = 1


f:id:ryooji_f:20170614214110p:image

ルジャンドル多項式が持つ性質の一つとして,”閉区間[−1, 1]上のL2-内積に関して直交する”というものがあります

%o10式に示す内積積分を 10 x 10 まで具体的に計算した結果が%o11式となります

これより%o10式が直交性を持つこと,対角成分は正規化されておらず 2/(2*n+1) となることが判ります

2017-04-27

佐藤・テイト予想

以前に谷山・志村予想のお話をした流れで,今回は佐藤・テイト予想(Sato–Tate conjecture)について触れておきたいと思います

この予想もRichard Taylorらによって既に証明されています


ここでは,難波完爾氏のテキストにある例題をフォローしてみます

sato-tate conjecture.wxm

http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170427/20170427004138.png

デデキントのη関数(Dedekind Eta function)として%o2式を定義します

η関数を使ってq展開式を与えます(%o3)

これをq^100まで展開した結果を%o4式に示します


f:id:ryooji_f:20170603152911p:image

n : N/2

いま素数 p = 2*n+1 に対してq展開式のq^nの係数をaとすると,p = 101までのpとaの組み合わせは%t6〜30となります

ここで,2次方程式 x^2-a*x+p = 0 の解を配列Cに代入していきます(画面出力は省略)



http://cdn-ak.f.st-hatena.com/images/fotolife/r/ryooji_f/20170603/20170603152906.png

N = 10000(p = 9973)までの解(複素根)を複素平面上にプロットしたものが%t1235です

これより解の分布が偏角(argument)θに依存していることが伺えますが,これがsin2θに比例するというのが同予想になります


佐藤幹夫の数学 増補版

佐藤幹夫の数学 増補版