Hatena::ブログ(Diary)

Maximaでこうぞうりきがく

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の組み合わせは%t32〜56となります

ここで,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θに比例するというのが同予想になります


佐藤幹夫の数学 増補版

佐藤幹夫の数学 増補版

2017-03-24

群 その2

前回に引き続いて群のお話です

構造力学に関係のありそうな具体的な演算について,群の定義を満足するかどうかをmaximaで確認してみます

group2.wxm

f:id:ryooji_f:20170324210050p:image

ベクトルの平衡移動について考えます

平面上のベクトルaの成分を%o1式で定義します(3番目の成分の"1"はダミーです)

e1方向にΔ1,e2方向にΔ2だけ平衡移動させる線形変換Pを%o2式で定義します

P.a はベクトルなので閉じています(%o3)

結合法則が成り立ちます(%o4)

単位元単位行列 I)が存在します(%o5, 6)

逆元(P^-1)が存在します(%o7, 8)

ここで,P^-1はPの反対方向への平衡移動に相当します


ということで,ベクトルの平衡移動は群を成します


f:id:ryooji_f:20170324210049p:image

ベクトルの鏡像反転について考えます

ベクトルaの成分を%o9式で定義します

e2-e3平面に対して鏡像反転させる線形変換Mを%o10式で定義します

M.a はベクトルなので閉じています(%o11)

結合法則が成り立ちます(%o12)

単位元単位行列 I)が存在します(%o13, 14)

逆元(M^-1)が存在します(%o15, 16)

ここで,M^-1はM自身であり,2回の操作で基に戻ることが分かります


ということで,ベクトルの鏡像反転は群を成します


f:id:ryooji_f:20170324210048p:image

ベクトルの回転変換について考えます

e3軸周りにθ[rad]だけ回転させる線形変換Rを%o17式で定義します

R.a はベクトルなので閉じています(%o18)

結合法則が成り立ちます(%o19)

単位元単位行列 I)が存在します(%o20, 21)

逆元(R^-1)が存在します(%o22, 23)

ここで,R^-1はRの逆周り-θ[rad]の回転変換に相当します


ということで,ベクトルの回転変換は群を成します