大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

仮想ばね法により基本境界条件を緩和させたB-spline Ritz法に設定する最適な仮想ばね定数の値について（その２）－自重の作用による均質な等断面長柱の線形座屈問題－

山本 寧音1・名木野 晴暢2†・足立 忠晴31長岡技術科学大学大学院，2都市・環境工学科，3豊橋技術科学大学†連絡著者 / Corresponding author (E-mail: nagino@oita-ct.ac.jp)

本稿では，線形化された梁－柱理論に従う自重の作用による均質かつ真っすぐな長柱の線形座屈問題を仮想ばね法により基本境界条件を緩和させたB-spline Ritz法により解析し，座屈荷重に与えるspline次数，解析領域内に配置するノットの数，仮想ばね定数の値および支持条件の影響を詳細に調べた．数値解に影響を与えない範囲の無次元仮想ばね定数の最大値は，前報と同様の鉛直ばねおよび回転ばね共に1010であった．この値を用いて，spline次数を5次に設定したときの座屈荷重は，支持条件にかかわらずノットの数の増
大にともなって有効桁数八桁から九桁で収束した．また，片持長柱の座屈荷重は有効桁数八桁から九桁まで厳密解に一致し，これ以外の支持条件については先行研究よりも精度の高い上界と思われる値を示した．

キーワード : Ritz法，正規化されたB-spline基底関数，仮想ばね法，
線形座屈問題，軸方向分布力，線形化された梁－柱理論

１．緒言ばその強形式は変数係数の四階線型の斉次常微分方程式の境界値問題で表される．これは，柱の断面諸量や材料特第一報1)では，線形化された梁－柱理論2)に従う中心軸圧性が部材軸方向に変化することによる．文献5)には，幾つ縮力を受ける均質かつ真っすぐな長柱（以下，長柱と略す）

かの変断面長柱と軸方向に不均質な長柱の線形座屈問題の線形座屈問題を仮想ばね法により基本境界条件を緩和させたB-spline Ritz法（以下，本手法と略す）により解析し，の厳密解が纏められている．これは第一種および第二種
Bessel関数により表現されるため，その数量化は容易でな

座屈荷重の精度，離散化条件，仮想ばね定数の値および支

く，ある程度の専門知識を必要とする．また，構造設計や持条件の関係を詳細に調べた．調査の結果，支持条件とspline次数にかかわらず，ノットの増大にともなって数値

材料設計の観点に立てば，柱の軸方向に任意に変化する断的に安定した近似解を与える仮想ばね定数の最大値は，鉛厳密解を求めることはほぼ不可能になり，何らかの手段に

10

面諸量や材料特性の取り扱いが求められる．こうなると，直ばねと回転ばね共に10 であった．また，関数空間に5次より強形式を近似的に解かざるを得ない．近似解には，対の正規化されたB-spline基底関数（以下，B-spline基底関数

象としている工学的諸問題を解決できるだけの精度が要と略す）を張ったときの本手法の座屈荷重はノットの数の求されることは言うまでもない．他方で，渡辺6)が述べて増大にともなって上から厳密解に近づき，収束値は有効桁いるように，何らかの手法により数値的に安定して高精度数八桁から九桁まで厳密解に一致するという高い精度をな近似解が得られるようになれば，それは厳密解や級数解示した．この問題の強形式は，定数係数の四階線型の斉次に代わって有限要素の性能評価や各種近似解析手法が強常微分方程式の境界値問題で表されるので，閉じた形の厳形式を正しく解いているかの検証へ活用できるであろう．

密解を求めることができ，その数量化は比較的容易である．

このような背景の下，本稿では変数係数の四階線型の斉これに対して，構造物の合理的な設計や材料の不均質性

次常微分方程式の斉次境界値問題で記述される“自重によの活用の見地から用いられる変断面長柱3)や軸方向傾斜機

る均質かつ真っすぐな長柱の線形座屈問題”を対象として，

能材料(Axially Functionally Graded Materials: AFGMs)4)の線

第一報1)と同様の調査を実施し，ある程度の有益な情報を

2)

形座屈問題は，線形化された梁－柱理論 に基づくとすれ得ることができたので報告する．

―38―

大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

２．数学モデル（強形式）

応力を軸線に関する一般化応力（断面力）として表したものである．式(1b)の下線部は，軸力N (x)が影響する非線形(1) 座標系と仮定項であることを意味している．この項を無視すると，式(1b)は初等梁理論のz軸方向の力のつり合い条件式に一致する．
長柱の左端点(x = 0)と右端点(x = L)での境界条件式は，

L
p (x)
x, ux

O

次のような組み合わせにより表される．

u = u on u or

nx N = N on ,(2a)w = w on w or

nxQ = Q on Q,(2b)on  ornx M = M on M.

(2c)

y, uy
z, uz

図-1 直線直交座標系O-xyzと長柱q = q

ここで， = {0, L}は一次元領域の境界を表しており，下付
図-1に示すような右手直交直線座標系O-xyzにある長柱き文字のu, w, は変位の境界条件を，N, Q, Mは一般化応力の単位長さ当たりの軸方向分布力p (x)の作用による分岐座の境界条件を意味している．nxは線部材の境界に立てた外屈問題を考える．ここで，x軸，y軸，およびz軸の設定は前

向き法線の方向余弦のx方向の成分であり，左端点(x = 0)の報1)と同様であり，ux (x, y, z), uy (x, y, z), uz (x, y, z)は，それ境界でnx = – 1, 右端点(x = L)の境界でnx = + 1を取る．はぞれ，長柱のx, y, z方向の変位である．柱の長さLはyz断面

曲げたわみ曲線wの接線と軸線とがなす角を表すたわみ角の寸法に比べて卓越しているとし，yz断面の断面積をA，中であり，これは微小角(| | << 1)であるとすれば，次のよう立軸に関する断面二次モーメントをIで表す（正しくはIyでに表すことができる．

ある）．長柱はx軸方向に一様断面とし，かつ同方向に均質q(x ) =な材料からなるとする．
本稿では，軸方向分布力p (x)が作用することにより長柱dw(x ).
dx(3)が僅かにy軸まわりに曲げ変形を生じてつり合った状態をまた，Qは z軸に平行なせん断力であり，次のように表さ考える．この不安定現象の数学モデルには，線形化されたれる．

梁－柱理論2)を用いる．この理論に設けられる仮定，およQ(x ) =

びその仮定に基づく変位場とひずみ場についても前報1)を参照されたい．

dM (x )
dw(x ).
+ N (x )
dx
dx

(4)

これは，式(1b)の{ }の中の式に相当する．僅かに曲がった
長柱の断面の接線方向に生じるせん断力S (x)とは異なるこ(2) 一般化応力の平衡方程式と境界条件式
単位長さ当たりの軸方向分布力p (x)を受ける長柱が僅かに曲げ変形を生じてつり合ったとき，長柱のある点の近傍
における局所的な力のつり合いを考える．
領域 = (0, L)にある長柱の任意の点の近傍の微小な線分要素におけるx軸方向とz軸方向の力のつり合い条件式，とに注意する．なお，各条件式の右辺のバーを付した変位
成分u, w，たわみ角および一般化応力N, Q, M は，与える既知の値であることを意味している．
長柱の支持条件の数学モデルは，これらの境界条件式を組み合わせて，次のように理想化される．
ピン支点（ヒンジ支点）：

およびxz平面内の任意の点に関する力のモーメントのつw = 0 on w,

り合い条件式は，それぞれ，次のように表される．

dN (x )+ p(x ) = 0 ,
dx

(1a)

d ìïï
dw(x ) üïïíS (x ) + N (x )
ý = 0,
dx ïîï
dx ïþï

(1b)dM (x )- S (x ) = 0 in W .
dx

(1c)

ここで，N (x)は僅かに曲がった長柱の断面の法線方向に生M = 0 on M.

(5)

w = 0 on w,

q = 0 on .

(6)

M = 0 on M,

Q = 0 on Q.

(7)

固定端：

自由端：

これらは，斉次の境界条件式である．
(3) 一般化応力と変位で表された構成方程式
長柱に対して，本稿のように右手直交直線座標系O-xyz

じる軸力，S (x)は同断面の接線方向に生じるせん断力であり，M (x)はy軸まわりの曲げモーメントである（正しくはを設定すると，yz断面に設定した図心軸まわりの断面一次Myである）．これらは，yz断面に分布する垂直応力とせん断モーメントは零になる．このことから，一般化応力と変位―39―

大分工業高等専門学校紀要

第 56 号(令和元年 11 月)で表された構成方程式に含まれる軸方向変位u (x)と面外変

自由端：位w (x)に関する項は非連成化され，微小ひずみ・微小回転に表される．

-EI2ü
ìïï d u (x ) 1 é d w (x ) ù ïï,ú ïN (x ) = EA ïí+ êêýï2 ë d x úû ïïïï d x
ï
î
þ

(8a)

d 3w (x )dw (x ) ,Q (x ) = -EI+ N (x )
dx
dx 3

(8b)

2

M (x ) = -EI

d w (x ) .
dx 2

本問題は，変数係数の四階線型の斉次常微分方程式のる．求められる固有値が座屈荷重qcrであり，これに対応する固有関数wが座屈したときの曲げ変形を表す．

(8c)３．強形式に等価な汎関数の停留値問題

係数である．なお，式(8a)の下線を付した右辺の第二項を無視したものが，線形化された梁－柱理論の構成則になる．

仮想ばね法により基本境界条件を緩和させた変分原理の汎関数（全ポテンシャルエネルギ）は，次のように表される．

(4) 支配方程式と境界条件による強形式
式(1c)を式(1b)に代入し，更に一般化応力と変位で表され

P = U +V

た構成方程式の式(8c)を代入すれば，変位で表された平衡+方程式である長柱の線形座屈問題の支配方程式（基礎方程+式）が次のように得られる．
2

(9)

aw 0
2
bq 0

(w - w 0 )2

x =0

2

(q - q0 )

2

+

+

x =0

awL
2

bq L

(w - w L )2

.

(15)

x =L

2

(q - qL )2

x =L

ここで，Uは長柱に蓄積されるひずみエネルギ，Vは軸方向の等分布力（自重）p (x) = p0による外力ポテンシャルであ

さて，x軸方向のつり合いの微分方程式(1a)は，線部材の単り，それぞれ，式(16)と式(17)によって与えられる．

純引張理論の軸力N (x)に関する平衡方程式を表しており，

LU =曲げとは非連成である．よって，境界条件式(2a)が与えられれば，N (x)の厳密解は容易に求められる．長柱の自重を

L

V =-

は密度であり，gは重力加速度である．これを式(1a)に代入
し，長柱の右端点(x = L)で自由端の条件(N = 0 on N)の下でうに求められる．

(16)0

p (x) = p0である．ただし，p0 = gAは分布力の荷重強度，

解けば，自重の作用により生じる長柱の軸力N (x)が次のよ

2

EI ó æç d2w ö÷÷
ô çç
dx .
2÷
÷
2 ô
õ çè dx ÷ø

単位長さ当たりに作用する軸方向の分布力p (x)で表すと，N (x ) = -p0 (L - x ) ( 0 £ x £ L ).

(14)

d3w
dw
- p0 (L - x )= 0 on Q.
dxdx 3

斉次境界値問題の自明でない解を求める固有値問題であ

ここで，EAは伸び剛性，EIは曲げ剛性であり，Eは縦弾性d éêd2w(x ) ùúd éêdw(x ) ùú-EI+= 0 in W .
N (x )
2 ê
2 ú
ê
dx úû
dx ëê
dx ûú dx ëd2w
= 0 on M ,
dx 2

-EI

の有限変位に基づく梁－柱理論の構成方程式は，次のよう

2æ dw ö
p0 ó
ô (L - x )ççç ÷÷÷ dx .
çè dx ø÷
2 ô
õ0

(17)

また，w0, wLは，それぞれ，長柱の左端点(x = 0)と右端点(x = L)の面外変位wに抵抗する仮想の鉛直ばね定数，0, Lは，それぞれ，長柱の左端点(x = 0)と右端点(x = L)のたわみ(10)角に抵抗する仮想の回転ばね定数，式中の右辺のバーを式(10)を式(9)に代入すれば，自重の作用を受ける均質かつ付した面外変位wとたわみ角は境界での既知の値であり，等断面長柱の線形座屈問題の支配方程式は，次のように表その下付き文字の0とLは，それぞれ，x = 0, x = Lでの条件

される．

であることを意味する．

ïï
d4w (x )d ìdw (x ) üïï-EIíp0 (L - x )ý = 0 in W .
4
ï
x
dx ï
d
dx
ï
ï
î
þ

式(12)と式(13)の基本境界条件は斉次であるので，本稿(11)で考えるべき汎関数は，次のように表される．
L

境界条件式を面外変位wで表すと，
ピン支点（ヒンジ支点）：

w = 0 on w,

2

dw
-EI
= 0 on M.
dx 2

+

(12)

+

固定端：

w = 0 on w,

dw= 0 on .
dx

L

2

2
æ dw ö÷
p0 ó
EI ó æç d2w ö÷÷
ç
ô
÷
ô
ç
P=
ç 2 ÷÷ dx - 2 ô (L - x )ççç dx ÷÷ dx
2 ô
è
ø
õ
õ çè dx ÷ø

(13)

aw 0
2
bq 0
2

0

w
q2

2
x =0

x =0

+

+

awL
2

b qL
2

0

w

q2

.

2

(18)

x =L

x =L

この汎関数が停留するような微分可能，かつ連続な関数wを見出せば，それは強形式の弱解を得ることに等しい．

―40―大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

４．仮想ばね法により基本境界条件を緩和させた
B-spline Ritz法による離散化された座屈方程
式と座屈条件式の定式化

nCP = m + p - 1 .

(22)式(19)を式(18)の汎関数に代入すると，は有限の自由度で近似される．この近似された汎関数の停留条件：

式(18)の汎関数の要請をすべて満足する微分可能，か¶P= 0 (i = 1,2, , nCP - 1, nCP ) ,¶Ai

つ連続な変関数（長柱の面外変位）w (x)を見出すことは困難を極める．よって，変関数wを次のような有限自由度の試行関数で近似し，汎関数が停留するような近似解を探(23)より，次の離散化された座屈方程式を得る．

求する．

[(K0 + KP ) - p0KG ]A = 0 .


nCP

w(x ) @ å AN(x ) ,
i i, p

(24)

K

(19)

i =1

ここで，

= Ni , p A

LNi, p = [N 1, p (x ), N 2, p (x ), , N n -1, p (x ), N n , p (x )] ,
CP

CP

2

T

2

ó d Ni, p d N j , p
K0 = EI ôdx ,õ0 dx 2 dx 2

ここで，
(20)

(25)

は自由境界を有する長柱の剛性行列，
A = {A1, A2 , , An -1, An }T ,
CP

(21)

CP

KP = aw 0 NiT, p N j , p

であり，Tは転置を意味する．Ni, p (x)は次数pのB-spline基底+bq 0

関数，iは基底関数の番号，Aiはi番目の基底の重み係数であは，(m + 2p)個である．ただし，mは領域内（境界を含+ awL NiT, p N j , p

dNiT, p dN j , pdx

dx

り，nCPは基底関数の数である．領域に1の分割(partition ofunity)を満たすB-spline基底を張るために必要なノットの数

x =0

+ bqL
x =0

x =L

,

dNiT, p dN j , p
dx

dx

(26)

x =L

は基本境界条件を緩和するための仮想ばねの剛性行列であり，K = K0 + KPは長柱の剛性行列である．また，む）に設定するノットの数であり，2pは領域の外側（境LTdNi , p dN j , pó
KG = p0 ô (L - x )
dx ,
dx
dx
õ0

界を含まない）に設ける付加ノットの数である．基底関数の数nCPと境界を含む領域内に設定するノットの数m(27)との間には，次の関係が成立する．

x, u

x, u

x, u

x, u

L

L

O

L

L

z, w

z, w
O
(a) P-P(b) C-Fz, w
O(c) C-C図-2 数値実験に用いる長柱の支持条件

―41―

z, w
O
(d) C-P

大分工業高等専門学校紀要

第 56 号(令和元年 11 月)は長柱の幾何剛性行列であり，0は零の列ベクトルである．

式によって求めることができる7)．

剛 性 行 列 K0 と 幾 何 剛 性 行 列 KG は ， (p + 1) 点 の Gauss-

L

Legendreの数値積分により求めた．
式(24)の自明でない解は，次式の座屈条件式：

det[(K0 + KP ) - p0KG ] = 0 ,

pcr =
(28)

2

æ d2w ö1ó
ô EI ççç 2 ÷÷÷ dx
çè dx ÷ø
2ô
õ
0
L

2

æ dw ö
p0 ó
ô (L - x )ççç ÷÷÷ dx
ôçè dx ÷ø
2 õ

.

(33)

0

から固有値p0を求め，それを式(24)に代入して固有値p0にここで，分子は座屈による曲げ変形によって長柱に貯えら対応する固有ベクトルAを決定することで得られる．

れるひずみエネルギ（内力仕事）であり，分母は長柱の曲げ座屈変形により生じる軸圧縮変位と軸方向等分布力とがなす外力仕事を意味する．式(33)の右辺の表示は，線形

５．数値実験および考察化された梁－柱理論に基づく長柱の線形座屈問題におけここでは，軸方向等分布力p (x) (= p0)が作用する均質かるRayleigh商である7)．これより，座屈荷重pcrは，外力仕事

つ等断面長柱の線形座屈問題を本手法により解析し，座屈に対するひずみエネルギの大きさとして評価されると解荷重pcrに与えるspline次数p，境界を含む領域内に設定する

釈することができよう．なお，許容である試行関数を仮定

ノットの数（以下，ノットの数と略す）m，仮想ばね定数

したRitz法による近似解は，厳密解の上界を与える8), 9)ので，の値および長柱の支持条件の影響を詳細に調べる．

次の関係が成り立つ．

本稿では長柱の座屈荷重pcr，仮想の鉛直ばね定数w0, wL

pcr ³ pcrExact .

および仮想の回転ばね定数0, Lを，次のように無次元化して表している．

(34)

ここで，不等式の左辺はRitz法による近似解であり，右辺

pcrL3

*

p =

EI
3

a* =

b* =

aw 0L
EI

bq 0L
EI

は厳密解を意味する．ただし，前報1)によれば，関数空間に(29).

ットの増大にともなって有効桁数九桁の収束値を得るが，数値計算上の誤差と思われる原因により，その値は厳密解

3

=

=

高次のB-spline基底を張ったとき，本手法の座屈荷重はノ

awLL

,

EI
bqLL
EI

(30)

の下界を与えることがある．
さて，前報1)では，四つの支持条件について全て閉じた形式の厳密解を求めることが可能なため，本手法の近似解

.

(31)

また，前報の知見1)を参考にして，数値実験に用いる無次元仮想ばね定数*を次のように定義した．
k* = a* = b * .

の上界と下界を明確にすることができた．本稿が対象としている問題についても厳密解は示されているが，前報のように閉じた形式で与えられるのはC-F（片持）の支持条件に限定されているようであり5), 10)，その他の支持条件につ(32)

いては境界条件を満足するべき級数解の形式（以下，べき級数解と略す）で与えられている11) – 13)．そのため，上界と数値実験に設定する支持条件は，図-2に示す四ケースと下界の明確な区別が難しい．そこで，本稿では，次の手順

し，例えばC-Fのように表す．これは柱の下端(x = 0)で固定

により，本手法の近似解に影響を与えない範囲の最大の無Cであり，上端(x = L)で自由Fであることを意味する．また，

次元仮想ばね定数*の値と，これを設定したときの本手法

ピン支点（ヒンジ支点）はPで表す．なお，(a)と(c)は対称の近似解の精度を調べることにした．
[1] 厳密解が得られるC-F（片持）について前報と同様のの支持条件，(b)と(d)は非対称の境界条件であり，(a)と(b)は静定，(c)と(d)は不静定である．

調査を行い，結果を整理する．

通常，工学分野の諸問題の解決を目的とするならば，強[2]形式の近似解は有効数字三桁から四桁程度の精度があれ

残り三ケースの支持条件については，前報と[1]の本手法の近似解の収束性と計算精度を判断材料とし，Huang・Dareing11)のべき級数解と杉山ら12)のべき級数ば十分であるが，本稿は前報と同様の観点1)の下，無次元
座屈荷重p*を有効桁数十桁程度で整理した．なお，数値計

解の収束性と有効桁数六桁の収束値と比較すること

算には64 bit OSのpersonal computerを用い，Fortranの倍精度で，本手法の近似解の収束性と精度を考察する．

で計算した．

なお，Huang・Dareing11)と杉山ら12)は，自重を受ける長柱の支配方程式の面外変位を式(35)のようなべき級数で近似し，(1) Timoshenkoのエネルギ法に基づく座屈荷重の評価

支配方程式と境界条件から未定係数bn に関する連立方程Timoshenokoのエネルギ法によれば，軸方向等分布力p (x)(= p0)の作用による均質かつ等断面長柱の座屈荷重pcrは，次式を導出している．これを数値的に解くことで，座屈荷重と，これに対応する座屈モードの数値を得ている．

―42―

大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

表-1 自重の作用による長柱の無次元座屈荷重pに与えるspline次数p, ノットの数mと仮想ばね定数の影響: C-F
p



4

108

m

1010

1012

108

5

1010

1012

Mode 1Mode 2Mode 3Mode 411

7.83734735

55.9771972

148.517011

285.588375

21

7.83734721

55.9770304

148.508401

285.451764

31

7.83734720

55.9770282

148.508302

285.450336

41

7.83734720

55.9770280

148.508295

285.450236

517.83734720

55.9770280

148.508294

285.450221

61

7.83734720

55.9770280

148.508294

285.450217

117.83734756

55.9771987

148.517015

285.588383

21

7.83734744

55.9770320

148.508406

285.451772

31

7.83734744

55.9770299

148.508307

285.450344

41

7.83734744

55.9770297

148.508299

285.450244

51

7.83734744

55.9770297

148.508298

285.450229

61

7.83734743

55.9770297

148.508298

285.450226

11

7.83734834

55.9772044

148.517031

285.588418

21

7.83734706

55.9770293

148.508398

285.451758

317.83734747

55.9770301

148.508307

285.450346

41

7.83734723

55.9770282

148.508296

285.450237

51

7.83734749

55.9770301

148.508299

285.450231

61

7.83734741

55.9770295

148.508298

285.450225

11

7.83734720

55.9770297

148.508592

285.463595

21

7.83734720

55.9770280

148.508294

285.450238

31

7.83734720

55.9770280

148.508294

285.450216

41

7.83734720

55.9770280

148.508294

285.450216

517.83734720

55.9770280

148.508294

285.450215

61

7.83734721

55.9770280

148.508294

285.450215

117.83734747

55.9770317

148.508597

285.463604

217.83734744

55.9770297

148.508299

285.450247

31

7.83734743

55.9770296

148.508298

285.450225

41

7.83734744

55.9770297

148.508298

285.450224

51

7.83734743

55.9770297

148.508298

285.450224

61

7.83734744

55.9770297

148.508298

285.450224

11

7.83734943

55.9770458

148.508634

285.463663

21

7.83734782

55.9770325

148.508307

285.450264

31

7.83734737

55.9770292

148.508297

285.450222

41

7.83734742

55.9770296

148.508298

285.450224

51

7.83734731

55.9770288

148.508296

285.450219

61

7.83734745

55.9770297

148.508298

285.450224

Exact5)

7.83734744

55.9770297

148.508298

285.450224

Huang and Dareing (1969)11)

7.837355.977148.51

Sugiyama et al. (1977)12)

7.83735

55.9770

The structural mechanics handbook (1986)15)

7.8373

Wang and Ang (1988)

14)

7.84

Duan and Wang (2008)13)12)

7.8373*

Note: 杉山ら のMode 1のp は下から収束し，N = 24としたときに有効桁数六桁の収束値を得ている．Mode 2は上から収束し，N = 36としたときに有効桁数六桁の収束値を得ている．

―43―

大分工業高等専門学校紀要

第 56 号

N(令和元年 11 月)

ただし，

w(x) = åbn x .
n

(35)

n =0

z=

ここで，は無次元座標の変数であるが，その定義は文献

32 q crL ,
3 EI

(37)

11) と 文 献 12) で 異 な る こ と に 注 意 さ れ た い ． Huang ・である．ここで，J– 1/3 (z)は(– 1/3)次の第一種Bessel関数であDareing11)が示した数値は，近似項数Nに対する座屈荷重のる．この方程式を解析的に解くことは困難であるため，こ収束性が不明であるが，一次から三次までの座屈荷重を有

れを有効桁数十五桁の精度で数値的に解くこととし，次の効桁数五桁で与えてくれている．杉山ら12)は，近似項数N四つの零根を得た．

に対する一次と二次の座屈荷重の変化の情報を有効桁数

z 1 = 1.86635085887389,

六桁で提供してくれている．この他にもDuan・Wang13)の変

z 2 = 4.98785323143515,

数rn/3 (n = 0, 1, 2, 3)と一般化された超幾何関数F (a1, a2; b1,z 3 = 8.12426538193969,

b2, b3; r)との積を基底とした境界条件を満足する厳密解（厳(38)z 4 = 11.2635148254276.

密解は四つの項で表され，一般化された超幾何関数の定義には無限級数が含まれる）も提案されているが，文献11)とここで，zの下付添字は，座屈モードの次数を意味する．式文献12)よりも数値解に関する情報が少ない．

(38)を式(37)に代入すれば，一次から四次までの無次元座
屈荷重p*の厳密解の数値を求めることができる．

(2) 片持(C-F)長柱表-1より，spline次数pにかかわらず，近似解に影響を与

表-1は，支持条件がC-Fである自重の作用による長柱のえない最も大きな無次元仮想ばね定数*の値は，前報1)と一次から四次までの無次元座屈荷重pに与えるspline次数p,同様に1010である．* = 108の近似解は柔らかく，* = 1012

ノット数mと無次元仮想ばね定数 の影響を表形式で整理のそれは数値的に安定していない．Spline次数p = 5に設定したものである．分岐座屈荷重に相当するのは最低次の座

すると，一次の無次元座屈荷重pは有効桁数八桁で，二次屈荷重であるが，二次以上の座屈荷重とこれに対応する座から四次までのそれは有効桁数九桁での収束値を得てい屈モードは，座屈荷重を高めるための重要な情報を提供しる．ただし，二次の無次元座屈荷重pは，ノットの数mを21てくれる1)ことから，本稿でも一次から四次までの座屈荷から41に変化させる過程で数値が振動しており，かつm =重について検討した．無次元仮想ばね定数*は，108, 109,31のときに下界である．上界定理から判断すると，九桁目1010, 1011, 1012の五ケースについて検討したが，本稿では108,の数値の信頼性は低いと考えられる．なお，一次と二次の1010, 1012の結果を提供する．B-spline基底関数の次数pは，無次元座屈荷重p は厳密解と有効桁数八桁まで一致して前報と同様の理由により，p = 4, 5に設定した1)．Spline次数

おり，三次と四次のそれは有効桁数九桁まで一致する．前p = 4のとき，面外変位wは三階導関数まで連続になるので，

報1)と比較すると，一次と二次の座屈荷重の精度が若干低

せん断力Q (x)は一次の区分的多項式で近似される．当然の下している．これは，長柱が自重によって座屈するときのことではあるが，spline次数p = 5に設定すると，面外変位w

曲げ変形の形が前報よりも複雑になっているためである

は四階導関数まで連続であるので，せん断力Q (x)の近似式と考えられる．

*は二次の区分的多項式になる．ノットの数mは，11から61まで変化させ，有効桁数九桁で表した無次元座屈荷重pの(3) 支持条件がP-P, C-CおよびC-Pである長柱値の変化を調べた．なお，表のboldの数値は近似解が収束

表-2から表-4は，それぞれ，支持条件をP-P, C-C, C-Pとしし，かつ厳密解と一致したことを意味する．また，近似解て表-1と同様の検討を行ったものである．ただし，本節では収束したが，その収束値が厳密解と一致しないものには， は厳密解を有効桁数九桁程度で数値化したものがないた
boldの数値に下線を付した．先にも述べたとおり，ここでめ，上界と下界を明確に区別することができず，また，近は厳密解と比較することで，近似解に影響を与えない無次似解の収束値が何に収束したのかを正確に判断すること

元仮想ばね定数 の最大の値の候補と，その精度を検証す*が難しい．よって，ここでは近似解が有効桁数九桁で収束る．表には厳密解の数値が記載されていれば十分であるが， した数値をboldで表すこととし，Huang・Dareing11)の一次か
参考までにHuang・Dareing11) と杉山ら 12) のべき級数解，
13)

14)

Duan・Wang の厳密解，Wang・Ang のRitz法による近似
15)

解，および構造力学公式集 の数値も併記した．

ら三次までの座屈荷重のべき級数解，および杉山ら12)のべき級数解の収束状態とその収束値と比較することで，本手法の近似解の精度を判断する．

厳密解の数値は，文献5)と文献10)を参考にして数値的に

表-2より，spline次数pにかかわらず，近似解に影響を与求めた．ここで，自重を受ける片持長柱の座屈方程式は，えない範囲の無次元仮想ばね定数*の最大値は1010である．

次式で与えられる5), 10)．

また，spline次数p = 5かつ無次元仮想ばね定数* = 1010に設J -1/3 (z ) = 0 .

(36)

定したときの本手法の近似解は，ノットの数mの増大にともない数値的に安定した収束状態を示しており，有効桁数―44―大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

表-2 自重の作用による長柱の無次元座屈荷重pに与えるspline次数p, ノットの数mと仮想ばね定数の影響: P-P
p



4

108

m

1010

1012

108

5

1010

1012

Mode 1Mode 2Mode 3Mode 411

18.5687256

86.4315979

196.318759

353.186380

2118.5687239

86.4308408

196.291038

352.751647

31

18.5687239

86.4308316

196.290779

352.748588

41

18.5687239

86.4308309

196.290760

352.748378

51

18.5687239

86.4308308

196.290757

352.748346

61

18.5687239

86.4308308

196.290757

352.748339

11

18.5687266

86.4316031

196.318775

353.186414

21

18.5687249

86.4308460

196.291055

352.751680

31

18.5687248

86.4308368

196.290796

352.748622

41

18.5687248

86.4308361

196.290777

352.748412

51

18.5687248

86.4308360

196.290774

352.748380

61

18.5687248

86.4308359

196.290773

352.748372

11

18.5687259

86.4316024

196.318774

353.186413

21

18.5687247

86.4308460

196.291055

352.751680

31

18.5687248

86.4308368

196.290796

352.748622

41

18.5687248

86.4308361

196.290777

352.748412

51

18.5687248

86.4308360

196.290774

352.748380

61

18.5687248

86.4308360

196.290773

352.748373

11

18.5687239

86.4308564

196.292612

352.791084

21

18.5687239

86.4308308

196.290761

352.748427

31

18.5687239

86.4308308

196.290756

352.748338

41

18.5687239

86.4308308

196.290756

352.748335

51

18.5687239

86.4308307

196.290756

352.748335

61

18.5687239

86.4308308

196.290756

352.748335

11

18.5687248

86.4308615

196.292629

352.791118

2118.5687248

86.4308360

196.290777

352.748460

31

18.5687248

86.4308359

196.290773

352.748371

41

18.5687248

86.4308359

196.290773

352.748369

51

18.5687248

86.4308359

196.290773

352.748369

61

18.5687248

86.4308359

196.290773

352.748369

11

18.5687245

86.4308614

196.292629

352.791118

2118.5687250

86.4308362

196.290777

352.748461

3118.5687248

86.4308360

196.290773

352.748372

41

18.5687248

86.4308360

196.290773

352.748369

51

18.5687249

86.4308360

196.290773

352.748369
352.748369

61

18.5687248

86.4308360

196.290773

18.56986.431

196.29Sugiyama et al. (1977)12)

18.5687

86.4309

The structural mechanics handbook (1986)15)

18.569Wang and Ang (1988)14)

18.58

Huang and Dareing (1969)

Duan and Wang (2008)
12)

11)

13)

18.5687*

Note: 杉山ら のMode 1のp は上から収束し，N = 30としたときに有効桁数六桁の収束値を得ている．Mode 2は下から収束し，N = 42としたときに有効桁数六桁の収束値を得ている．

―45―

大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

表-3 自重の作用による長柱の無次元座屈荷重pに与えるspline次数p, ノットの数mと仮想ばね定数の影響: C-C
p



4

108

1010

1012

108

5

1010

1012

m

Mode 1Mode 2Mode 3Mode 411

74.6290737

157.044090

325.758330

495.467677

21

74.6285668

157.032875

325.516057

494.004251

31

74.6285605

157.032742

325.513599

493.989459

41

74.6285601

157.032732

325.513431

493.988483

51

74.6285600

157.032731

325.513405

493.988339

61

74.6285600

157.032731

325.513400

493.988306

1174.6290822

157.044139

325.758384

495.467860

21

74.6285755

157.032924

325.516112

494.004432

31

74.6285692

157.032792

325.513654

493.989639

41

74.6285687

157.032782

325.513486

493.988663

51

74.6285687

157.032780

325.513460

493.988519

61

74.6285686

157.032780

325.513454

493.988487

1174.6290932

157.044167

325.758436

495.468037

21

74.6285710

157.032915

325.516094

494.004406

3174.6285705

157.032794

325.513659

493.989647

41

74.6285668

157.032778

325.513478

493.988652

51

74.6285694

157.032782

325.513463

493.988525

61

74.6285684

157.032780

325.513454

493.988486

1174.6285675

157.033185

325.542366

494.427873

2174.6285600

157.032731

325.513442

493.988741

31

74.6285600

157.032730

325.513398

493.988303

4174.6285600

157.032730

325.513397

493.988291

51

74.6285600

157.032730

325.513397

493.988290

61

74.6285600

157.032730

325.513397

493.988290

1174.6285766

157.033235

325.542422

494.428056

2174.6285687

157.032781

325.513498

493.988923

31

74.6285686

157.032780

325.513453

493.988483

41

74.6285686

157.032780

325.513452

493.988472

51

74.6285686

157.032780

325.513452

493.988471

61

74.6285686

157.032780

325.513452

493.988471

1174.6286030

157.033287

325.542508

494.428094

21

74.6285738

157.032792

325.513522

493.988964

31

74.6285671

157.032777

325.513447

493.988475

41

74.6285683

157.032779

325.513451

493.988471

51

74.6285673

157.032777

325.513446

493.988463
493.988473

61

74.6285688

157.032780

325.513453

Huang and Dareing (1969)11)

74.629157.03325.51

Sugiyama et al. (1977)12)

74.6286

157.031

The structural mechanics handbook (1986)15)

74.629Wang and Ang (1988)14)78.96

Duan and Wang (2008)
12)

13)

74.6286*

Note: 杉山ら のMode 1のp は上から収束し，N = 36としたときに有効桁数六桁の収束値を得ている．Mode 2は下から収束し，N = 42としたときに有効桁数六桁の収束値を得ている．

―46―

大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

表-4 自重の作用による長柱の無次元座屈荷重pに与えるspline次数p, ノットの数mと仮想ばね定数の影響: C-P
p



4

108

m

1010

1012

108

5

1010

Mode 1Mode 2Mode 3Mode 411

52.5007890

129.059252

276.364762

439.725890

2152.5006621

129.054333

276.245006

438.864630

3152.5006604

129.054274

276.243769

438.855769

41

52.5006603

129.054270

276.243682

438.855176

51

52.5006603

129.054269

276.243669

438.855088

6152.5006603

129.054269

276.243666

438.855068

1152.5007917

129.059270

276.364784

439.725970

21

52.5006648

129.054352

276.245029

438.864709

3152.5006632

129.054293

276.243792

438.855848

41

52.5006631

129.054288

276.243705

438.855255

51

52.5006631

129.054288

276.243692

438.855167

6152.5006630

129.054287

276.243689

438.855147

1152.5007961

129.059284

276.364814

439.726070

2152.5006626

129.054346

276.245016

438.864690

3152.5006633

129.054293

276.243793

438.855851

41

52.5006618

129.054285

276.243698

438.855245

51

52.5006634

129.054289

276.243694

438.855171

61

52.5006630

129.054287

276.243689

438.855147

1152.5006615

129.054416

276.254325

439.061301

2152.5006603

129.054269

276.243684

438.855290

31

52.5006603

129.054269

276.243665

438.855065

4152.5006603

129.054269

276.243665

438.855059

51

52.5006603

129.054269

276.243665

438.855058

6152.5006603

129.054269

276.243665

438.855058

1152.5006645

129.054435

276.254349

439.061382

2152.5006631

129.054288

276.243707

438.855370

31

52.5006630

129.054287

276.243688

438.855144

41

52.5006630

129.054287

276.243688

438.855138

51

52.5006630

129.054287

276.243688

438.855138

6152.5006630

129.054287

276.243688

438.855137

1152.5006776

129.054467

276.254408

439.061413

21

52.5006656

129.054295

276.243723

438.855399

3152.5006627

129.054286

276.243686

438.855141

41

52.5006630

129.054287

276.243687

438.855138

51

52.5006622

129.054286

276.243683

438.855131

61

52.5006631

129.054288

276.243688

438.855139

52.501129.05276.24

Sugiyama et al. (1977)12)52.5009129.054

The structural mechanics handbook (1986)15)

52.501

Wang and Ang (1988)14)

53.91

1012

Huang and Dareing (1969)11)

Duan and Wang (2008)
12)

13)

52.5007*

Note: 杉山ら のMode 1のp は上から収束し，N = 30としたときに有効桁数六桁の収束値を得ている．Mode 2は下から収束し，N = 42としたときに有効桁数六桁の収束値を得ている．

―47―

大分工業高等専門学校紀要

第 56 号(令和元年 11 月)いると考えられる．

九桁での収束値を得る．また，その値は先行研究とよく一1)致している．前報 と本稿の片持長柱の数値実験結果から

なお，固有関数とその導関数の精度についても検討してい判断すると，本手法の近似解は，先行研究よりも精度の高るが，これについては別の機会に報告する予定である．
前報 1) と本稿の調査から汎関数の変関数をB-spline基底

い上界の値を与えていると考えられるが，二次の座屈荷重についてはもう少し詳細な検討が必要である．

関数により近似するとき，仮想ばね法により基本境界条件

表-3からは，次のことがわかる．

を緩和させることによっても有効桁数八桁までは厳密解 Spline次数pにかかわらず，数値的に安定し，近似解にに一致する，または一致していると思われる高い解析精度影響を与えない最大の無次元仮想ばね定数*の値は，が得られることを示した．今後は，基本境界条件を正確に満足する許容な試行関数を仮定したB-spline Ritz法により10

やはり10 である．
 C-FおよびP-Pの支持条件に比べると，近似解の収束性
1)が良くない．これは，前報 と同様である．

同様の調査を実施し，弱形式の近似解にB-spline基底関数を用いたときの基礎的な情報を提供したいと考えている．

 Spline次数p = 5かつ無次元仮想ばね定数* = 1010に設定

なお，本稿は同研究分野の専門家による査読を経ずに公開したときの本手法の近似解は，ノットの数mの増大にとされるが，本稿で示した表の数値は，読者の判断で種々のもない数値的に安定した収束状態を示しており，有効離散化手法の検証などに使っていただければ幸いである．

桁数九桁での収束値を得る．その値は先行研究とよく一致している．
 P-Pの支持条件と同様，本手法の近似解は，先行研究よ

謝辞：本研究の一部は，LIXIL住生活財団調査研究助成（2013年度）および前田記念工学振興財団研究助成（平成

りも精度の高い上界の値を与えていると考えられるが， 27年度）を受けて行われました．ここに記して謝意を示し
ます．

二次の座屈荷重については検討が必要である．
表-4についてもこれまでとほぼ同様の結果が得られている．異なるのは，次の二点であった．
 Spline次数p = 5かつ無次元仮想ばね定数* = 1010に設定

参考文献
1)

したときの本手法の一次から三次の無次元座屈荷重p

名木野晴暢，山本寧音，足立忠晴：仮想ばね法により基本境界条件を緩和させたB-spline Ritz法に設定する

はノットの数mの増大にともなって有効桁数九桁で収最適な仮想ばね定数の値について（その１）－中心軸束するが，四次のそれは有効桁数八桁である．

圧縮力を受ける均質な等断面長柱の線形座屈問題－,
大分工業高等専門学校紀要, 第56号, pp.24-37, 2019.

 本手法の近似解は，先行研究よりも精度の高い上界の値を与えていると考えられるが，一次の座屈荷重につ

2)

６．結

土木学会構造工学委員会計算力学とその応用に関する研究小委員会：構造工学における計算力学の基礎と応用 (構造工学シリーズ 7), 土木学会, pp.118-154,
1997.

いてはもう少し詳細な検討が必要である．

言

本稿では，線形化された梁－柱理論に従う自重の作用に3)

崎山毅，栗林和夫：変断面梁柱の挫屈解法について，
長崎大学工学部研究報告, 第8号, pp.25-31, 1977.

4)

山本寧音，名木野晴暢，足立忠晴，樋口理宏：軸方向よる均質かつ真っすぐな長柱の線形座屈問題を仮想ばね法により基本境界条件を緩和させたB-spline Ritz法により傾斜機能材料からなる柱部材の線形座屈特性に与え解析し，座屈荷重の精度，spline次数，解析領域内に配置す学研究会シンポジウム論文集 (USB), R4-1 (6 pages),2017.

る縦弾性係数の変化の影響，第5回九州橋梁・構造工

るノットの数，仮想ばね定数の値および支持条件の関係を詳細に調べた．本調査から，次の知見が得られた．
[1] 支持条件とspline次数pにかかわらず，近似解に影響を5)buckling of structural members, CRC Press, pp.9-75(chapter 2), 2004.

与えない範囲の仮想ばね定数の最大値は，無次元仮想鉛直ばね定数*，無次元仮想回転ばね定数*共に

6)

10

Spline次数p = 5かつ無次元仮想ばね定数* = 1010に設

7)

荷重p*は，支持条件にかかわらず，ノットの数mの増8)

T.R.トーカート (著)．岩本卓也 (翻訳)：構造力学とエネルギ原理 (理工学海外名著シリーズ32), ブレイン

大にともなって上から正解に近づくような収束状態を示し，有効桁数八桁から九桁での収束値を得る．
[2]で得られた収束値は，有効桁数八桁から九桁の精

福本秀士：構造物の座屈・安定解析 (新体系土木工学
9), 技報堂出版, pp.37-98, 1982.

定したときの長柱の一次から四次までの無次元座屈[3]

渡辺力：p-Ritz法における試行関数の改良に関する一考察，構造工学論文集A, 61A巻, pp.1-10, 2015.

10 である．
[2]

C.M. Wang, C.Y. Wang, J.N. Reddy: Exact solutions for

図書出版, pp.76-77, 1979.
9)度で，先行研究よりも精度の高い上界の値を与えて―48―

C.L.ディム，I.H.シャームス (共著)，砂川恵 (監訳)：
材料力学と変分法 (理工学海外名著シリーズ26), ブ

大分工業高等専門学校紀要

10)

第 56 号(令和元年 11 月)

レイン図書出版, pp.94-151, 1977.
S.P. Timoshenko and J.M. Gere: Theory of Elastic Stability

13)columns including self-weight, Journal of Engineering
Mechanics (ASCE), Vol. 134 (1), pp. 116-119, 2008.

(Dover Civil and Mechanical Engineering), Dover
14)

Publications (2nd ed. edition), pp.100-113, 2009.
11)

C.M. Wang and K.K. Ang: Buckling capacities of braced

T. Huang and D.W. Dareing: Buckling and Frequencies ofheavy columns under an axial load, Computers &Long Vertical Pipes, Journal of the Engineering Mechanics

Structures, Vol. 28 (5), pp.563-571, 1988.

Division (ASCE), Vol. 95(1), pp. 167-182, 1969.
12)

W.H. Duan and C.M. Wang: Exact solution for buckling of

15)

杉山吉彦，芦田幸逸，川越治郎：自重による長柱の座
屈, 日本機械学會論文集, 43 巻, 376号, pp.4435-4443,
1977.

―49―

土木学会構造工学委員会構造力学公式集改訂委員会
編：構造力学公式集 (昭和61年版), 土木学会, p.118,
1986.
(2019.9.30 受付)