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

第 56 号(令和元年 11 月)

仮想ばね法により基本境界条件を緩和させたB-spline Ritz法に設定する最適な仮想ばね定数の値について（その１）－中心軸圧縮力を受ける均質な等断面長柱の線形座屈問題－

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

本稿では，線形化された梁－柱理論に従う中心軸圧縮力を受ける均質かつ真っすぐな長柱の線形座屈問
題を仮想ばね法により基本境界条件を緩和させたB-spline Ritz法により解析し，分岐座屈荷重の精度，spline
次数，解析領域内に配置するノットの数，仮想ばね定数の値および支持条件の関係を詳細に調査した．数
値的に安定し，近似解に影響を与えない範囲の無次元仮想ばね定数の最大値は，鉛直ばねおよび回転ばね
共に1010であった．この仮想ばね定数値を用いて，spline次数を5次，ノットの数を41に設定したときの本手法の座屈荷重は，支持条件にかかわらず有効桁数九桁で収束した値が得られた．また，その収束値は有効
桁数八桁から九桁まで本問題の強形式の厳密解に一致するという高い解析精度を示した．

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

１．緒言な物体の基本振動数の近似値を算定する方法）を，1908年にWalther Ritzが一般化した変分原理に基づく近似解析手

本稿は，線形化された梁－柱理論1)に従う中心軸圧縮力法である．Ritz法は，対象としている問題の強形式の正解を受ける均質かつ真っすぐな長柱の線形座屈問題を仮想（厳密解）を探す代わりに強形式をこれと等価な汎関数のばね法により基本境界条件（変位の境界条件またはDirichlet条件）を緩和させたB-spline Ritz法（以下，本手法

停留値問題に置き換え，この問題の変関数の正解を幾つかと略す）により解析し，座屈荷重の精度，離散化条件，仮を汎関数の停留条件から決定することで変分問題の近似

想ばね定数の値および支持条件の四者の関係を調べ，その解を得る離散化手法である．Ritz解は強形式の厳密解に要情報を提供するものである．本調査は，著者らが着手してい る 軸 方 向 傾 斜 機 能 材 料 (Axially Functionally Graded

請される事項の幾つかを満足しない弱解である．Ritz法にの既知関数の線形結合によって近似し，各関数の重み係数は，次のような特色がある7)．

2) – 5)Materials: AFGMs)の力学的特性に関する研究において，  定式化が容易である．

変数係数を有する四階線型の斉次常微分方程式の境界値

 汎関数の変関数の近似に許容な試行関数を構成できれ問題（強形式）を精度良く解くための離散化手法を確立さばよく，自然境界条件の処理は自動的に，無造作に実せる過程3)で実施した．この研究の背景や目的などについ行される．なお，この条件は近似的に満足される．
 強形式を汎関数の停留値問題に表現し直せるとき，Ritzては，文献2)を参照されたい．

法により得られる近似解はある意味で厳密解の最良近似となる．

２．B-spline Ritz法

 一般的な問題に適用できる重み付き残差法に対して，
Ritz法の適用は対象としている問題の強形式を汎関数(1) Ritz法の停留値問題に置き換えられるときに限定される．

Rayleigh-Ritz法6)（以下，Ritz法と略す）は，1877年にLord

 長方形領域などの問題には適用しやすいが，境界の形

Rayleighが Theory of Sound で提唱したRayleigh法（連続的状が複雑になると，基本境界条件を満足する近似関数―24―大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

ただし，Weierstrassの定理およびDu Bois-Reymondsの定理の構成が困難になる．
なお，強形式の厳密解や汎関数の停留関数の自由度は，一によれば，一般に汎関数に高階の導関数が含まれていても，般に無限大である．これに対してRitz法は，無限の自由度

導関数の区分的な連続性だけを仮定すれば，許容される関を有する正解を有限の自由度で近似する．この“自由度を数の条件を満たす9)．試行関数は，一般に線形独立な既知有限個に減じること”は，幾何学的拘束を余分に加えるこの基底関数（座標関数）i(x)の線形結合により表されることに相当する8)．

とが多い．特に，一次元問題を対象とすれば，汎関数のRitz法では，汎関数の変関数を“許容である試行関数”

変関数を近似する許容な試行関数f(x)は，次のように仮定される8), 10), 11)．

により近似する．ここで，許容である関数とは，次の二つの条件を満足する関数のことをいう．
[1]nf (x ) = f0 (x ) + åaifi (x ) .

近似関数は，汎関数に要求される最も高次の導関数
まで微分可能，かつ連続でなければならない．
近似関数は，基本境界条件を満足していなければな

ここで，ai (i = 1, 2, …, n)は，各基底関数の影響度を表す重らない．

み係数（展開係数）である．0(x)は，境界での基本境界条

Ni,1(t)

1.0

i=1

i=2

i=4

i=3

i=5

0.5

0.0
-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

1.25

1.50

1.75

2.00

t
(a) 一次の正規化されたB-spline基底関数Ni, 1 (t); p = 1, min = 5, nCP = 5, t = {- 1 , 0, 1 , 1 , 3 ,1, 5 }
4

4

2

4

4

1.0

Ni,2(t)

i=1

i=2

i=3

i=4

i=5

i=6

0.5

0.0
-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50
t

0.75

1.00

1.25

1.50

1.75

2.00

(b) 二次の正規化されたB-spline基底関数Ni, 2 (t); p = 2, min = 5, nCP = 6, t = {- 1 , - 1 , 0, 1 , 1 , 3 ,1, 5 , 3 }
2

4

4

2

4

4 2

1.0

Ni,3(t)

[2]

(1)

i =1

i=1

i=2

i=3

i=4

i=5

i=6

i=7

-0.25

0.00

0.25

0.50
t

0.75

1.00

1.25

0.5

0.0
-1.00

-0.75

-0.50

1.50

1.75

2.00

(c) 三次の正規化されたB-spline基底関数Ni, 3 (t); p = 3, m = 5, nCP = 7, t = {- 3 , - 1 , - 1 , 0, 1 , 1 , 3 ,1, 5 , 3 , 7 }4

図-1 次数pの正規化されたB-spline基底関数Ni, p (t)

―25―

2

4

4

2

4

4 2 4

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

第 56 号(令和元年 11 月)

件を満足するように適当に設定する．他方，n個の基底関数

基本境界条件を正確に満足させることができる．しかし，

i(x) (i = 1, 2, …, n)は，境界で斉次の基本境界条件を満足し，基本境界条件の数，すなわち，Lagrangeの未定乗数の数だかつ汎関数の中に現れる導関数の最高階数まで微分可能な連続関数を選ぶ．このとき，近似された系の汎関数 (f)け未知量が増大する．また，係数行列の対角要素に零が入は展開係数aiの関数でもあるので，汎関数の停留条件は，きない11)．後者は，前者と異なり未知量の数は変化せず，¶P
dai = 0 .
i =1 ¶ai

係数行列の正定値行列も保たれる．また，数学的には，仮ndP(a1, a2, , an -1, an ) = å

るので正定値行列を前提としたソルバーを使うことがで(2)

想ばね定数の値（penalty数）を無限大にすると変位は拘束され，その値を零にすれば自由表面となる12), 13), 17) – 19)．しここで，aiは展開係数aiの変分である．これは任意である

かし，数値計算においては無限大を数値として取り扱うこことから，次のn元連立代数方程式を得る．

とができないため，解に影響を与えない範囲の最大値を設定することになる．その為，数値解析に用いる仮想ばね定¶P= 0 (i = 1,2, , n ) .
¶ai

(3)

よって，式(3)のn元連立代数方程式を解けば展開係数aiが求められ，式(1)の近似関数が決定する．なお，仮定する近似関数は許容されるものであり，かつn個の基底関数列 i(x)が完備であれば，n → ∞ としたときの極限において，試行関数f (x)はエネルギの意味において考えている問題の正解に収束する10)．

数の値は，事前に，試行錯誤的に求めておかなければならない．文献12)と文献13)は，仮想ばね法により斉次の基本
境界条件を近似的に満足させているが，精度の高い近似解が得られている．
(3) Ritz法，B-spline Ritz法と有限要素法の関係
有限要素法は，偏微分方程式の初期値境界値問題の近似解を得るための強力な離散化手法の一つである．これとRitz法およびB-spline Ritz法との関係を簡単に述べておく．
一般に，領域の全体に既知の基底（全体基底関数）を(2) B-spline Ritz法
Ritz法の近似解の精度は，式(1)の0 (x)とn個の基底関数

i (x)の選択が重要な鍵となる．B-spline Ritz法12), 13)は，弾性問題の強形式を汎関数の停留値問題に変換したときの変関数（未知の変位）に正規化されたB-spline基底関数（以下，
正規化されたを省略する）と重み係数の積の線形結合で近似したRitz法である．Spline関数は，多項式を何らかの連続条件を満たすように接続した区分的な多項式である．Bspline基底関数は，spline関数を局所的な台を持つ基底によって表現したものである．文献11)と文献12)では，単純ノットベクトルtからなるB-spline基底関数を近似関数に選択している．パラメータ領域[0, 1]の内部に設定するノットの数mを5に固定し，単純ノットベクトルtからなる次数pのBspline基底関数Ni, p (t)を図-1に示した．ここで，tはパラメー
タ，tiは非減少の実数の列で定義されるノットであり，nCP
は基底関数の数である．領域[0, 1]の左右の外側にあるノッ
トは付加ノットであり，これに該当するノットベクトルtの張ったものをRitz法と呼び，領域を幾つかの小領域i (i =1, 2, …, n)の集合とし，その各小領域iに比較的低次の基底（局所基底関数）を張ったものを有限要素法と呼んでいる．
Ritz法は，領域を一つの要素とした有限要素法と解釈することもでき，有限要素法はRitz法の基底関数の選択を工夫したものに相当する．B-spline Ritz法は局所基底を張るので，分類としては有限要素法に近く，一次のB-spline基底を張ったときは，一次の有限要素に一致する．しかし，Bspline基底関数の有する高次の導関数までの連続性から，二次以上のB-spline基底を領域に張ると，一次のB-spline
基底と異なり必ずしも領域を小領域（要素）に“分割”iしているとは言い難い．よって，B-spline Ritz法は，Ritz法と有限要素法の両方の特徴を有する離散化手法であると考えられる．また，有限の多項式によるパラメータ表示によって物体の形状を正確に表現できないものを除けば，
isogeometric写像20), 21)も可能であると思われる．

成分には下線を付した．Ritz法の近似関数は，基本境界条
件を満足するように選ばなければならないが，図-1の基底関数は次数pに関係なく境界(t = 0とt = 1)で斉次の基本境界条件を満足しないことがわかる．ただし，図-1 (a)の基本境
界条件の処理は容易である．他方，図-1 (b)と(c)は，パラメ

３．数学モデル（強形式）(1) 座標系と仮定

ータtの境界(t = 0とt = 1)で有限の値を持つ基底関数がp個Lだけ存在するため，基本境界条件の処理は容易でない．こOのような近似関数を試行関数に選択するときは，Lagrangeの未定乗数法14) – 16)により基本境界条件を付帯条件とした変分原理で表現するか，または仮想ばね法 12), 13), 17) – 19)（penalty関数法）により基本境界条件を緩和させた変分原理で表現するか，に大別される．前者によれば，付帯する

―26―

y, uy
z, uz

図-2 直線直交座標系O-xyzと細長い線部材x, ux大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

図-2に示すような右手直交直線座標系O-xyzにある真っ直ぐな細長い線部材（以下，長柱と略す）の中心軸圧縮力
による分岐座屈問題を考える．図-2のux (x, y, z), uy (x, y, z),た，(dw/dx)は曲げたわみ曲線w (x)の接線と軸線とがなす角を表すたわみ角であり，微小角(| | << 1)とすることで，次式のように表すことができる．

uz (x, y, z)は，それぞれ，長柱のx, y, z方向の変位を意味する．
q(x ) =

Lは柱の長さであり，これはyz断面の寸法の長さに比べて

d w (x ) .
dx(5)

卓越しているとする．また，長柱はx軸方向に一様断面とし，
梁－柱理論のひずみ－変位の関係式は，次のように与えかつ同方向に均質な材料からなるとする．y軸とz軸は，長柱のyz断面の図心を通るように設定する．これにより，y軸られる．

は中立軸となる．x軸は，yz断面の図心を通るように設定すexx (x, y, z ) =をA，中立軸に関する断面二次モーメントをIで表す（正しくはIyである）．
本稿では，長柱の端部(x = L)に中心軸圧縮力が作用し，¶ux (x, y, z ).
¶x= e0 (x ) + z k0 (x )る（これを軸線と呼ぶ）．断面諸量として，yz断面の断面積

(6)

ここで，
2

これにより端部がx軸方向にuLだけ圧縮変位してつり合っ

e0 (x ) =ている系から，軸圧縮力の増加にともない長柱が僅かにydu(x ) 1 éê dw (x ) ùú ,+ êdx
2 ë dx úû(7)d2w (x )dq(x ) ,
=dx
d2x

(8)

軸まわりに曲げ変形を生じてつり合った系を考える．この不安定現象の数学モデルには，線形化された梁－柱理論1)k0 (x ) = -を用いる．これは，以下の仮定により構築される．
[1]

長柱の変位は，微小ひずみ，かつ微小回転とし，線形化された有限変位理論に従う．

[2]

変形前に軸線（x軸）に垂直であった断面は，曲げ座
屈変形後も軸線に垂直であるとする（Bernoulli-Euler
の仮定）．

[3]

であり， (x)は軸線方向の伸縮ひずみであり， (x)は僅かに曲がった長柱の断面の曲率を表している．
(3) 一般化応力の平衡方程式と境界条件式
中心軸圧縮力Pを受ける長柱が僅かに曲げ変形を生じて曲げ座屈変形後の断面形状不変を仮定する．

つり合ったとき，長柱のある点の近傍における局所的な力[4] 材料の構成則は，一軸応力状態のHookeの法則に従う．
のつり合いを考える．領域 = (0, L)にある長柱の任意の点
ここで，[2]の仮定は断面のゆがみの原因であるせん断変形の近傍の微小な線分要素におけるx軸方向とz軸方向の力を考慮しないためのものであり，xz面の面外せん断ひずみのつり合い条件式，およびxz平面内の任意の点に関する力
xz = 0，かつxy面の面内せん断ひずみxy = 0 であることを意のモーメントのつり合い条件式は，それぞれ，次のように味している．また，[3]の仮定は，y軸方向の垂直ひずみyy =表される．
0，z軸方向の垂直ひずみzz = 0，および yz面の面外せん断
dN (x )
ひずみyz = 0 であることに相当する．この仮定に従う長柱
(9a)
=0,
dx
の変形は，垂直ひずみによる長柱の伸縮の影響，変位による長柱の回転の有限性，および曲げ変形による部材軸力の変化の影響を無視した範囲に限定される1)．
(2) 変位場とひずみ場(1)の座標系と仮定に基づくとき，運動学的に許容される
長柱の任意の点(x, y, z)の変位ux (x, y, z), uy (x, y, z), uz (x, y, z)dw(x ),dx(9b)

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

(9c)

ここで，N (x)は僅かに曲がった長柱の断面の法線方向に生は，それぞれ，次式のように表すことができる．

ux (x, y, z ) = u(x ) - zd ìïïdw(x ) üïïíS (x ) + N (x )
ý = 0,
dx ïîï
dx ïþï

じる軸力，S (x)は同断面の接線方向に生じるせん断力であ(4a)り，M (x)はy軸まわりの曲げモーメントである（正しくはMyである）．これらは，yz断面に分布する垂直応力とせん断

uy (x, y, z ) = 0 ,(4b)uz (x , y, z ) = w(x ) .

(4c)

応力を軸線に関する一般化応力（断面力）として表したものである．式(9b)の下線部の項は，軸力N (x)が影響する非線形項である．この項を無視すると，式(9b)は初等梁理論のz軸方向の力のつり合い条件式に一致する．軸力N (x)がここで，u (x), w (x) は，それぞれ，断面の軸線(x, 0, 0)にお変数xの関数として表されるとき，式(9b)は変数係数の線型

けるx, z方向の並進変位であり，これらは独立した未知関数の斉次常微分方程式になり，厳密解の探求を複雑にするこである．以降，u (x)を軸変位，w (x)を面外変位と呼ぶ．まとは言うまでもない．

―27―

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

第 56 号(令和元年 11 月)

長柱の左端点(x = 0)と右端点(x = L)での支持条件，すな
Q (x ) = -EI

わち，境界条件には，一般化応力の境界条件式と変位の境界条件式がある．これらは，次のような組み合わせにより表される．

d 3w (x )d w (x ) ,+ N (x )
dx
dx 3

(15b)

d2w (x ) .
dx 2

(15c)

M (x ) = -EIu = u on u or nxN = N on ,(10a)

ここで，EAは伸び剛性，EIは曲げ剛性であり，Eは縦弾性

w = w on w or nxQ = Q on Q,(10b)nxM = M on M.

(10c)

q=q

on  or係数である．なお，式(15a)の下線を付した右辺の第二項を無視したものが，線形化された梁－柱理論の構成則になる．
(5) 支配方程式と境界条件式による強形式

ここで， = ∂ = {0, L}は一次元領域の境界を表しており，式(9c)を式(9b)に代入した後，更に一般化応力と変位で表下付き文字のu, w, は変位の境界条件（基本境界条件）を，された構成方程式の式(15c)を代入すれば，変位で表された

N, Q, Mは一般化応力の境界条件（自然境界条件）を表して平衡方程式である長柱の線形座屈問題の支配方程式（基礎いる．nxは長柱の境界に立てた外向き法線の方向余弦のx方方程式）が次のように得られる．

向の成分であり，左端点(x = 0)の境界でnx = – 1, 右端点(x =d2 êéd2w(x ) úùd éêdw(x ) ùú-EI+= 0 in W .
N (x )
2 ê
2 ú
ê
dx úû
dx êë
dx úû dx ëL)の境界でnx = + 1 を取る．また，Qは z軸に平行なせん断力であり，次のように表される．

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

(16)

さて，x軸方向のつり合いの微分方程式(9a)は，線部材の単(11)純引張理論（棒の理論）の軸力N (x)に関する平衡方程式に等しく，曲げに関する支配方程式式(16)とは非連成である．

これは式(9b)の{ }の中の式に相当し，僅かに曲がった長柱

よって，境界条件式(10a)が与えられれば，N (x)の厳密解はの断面の接線方向に生じるせん断力S (x)とは異なることに容易に求めることができる．式(9a)を長柱の右端点(x = L)で注意する．なお，各条件式の右辺のバーを付した変位u, w，

一般化応力の境界条件式(nxN = – P on N)で解けば，中心軸

たわみ角および一般化応力N, Q, Mは，境界条件として与圧縮力の作用により生じる長柱の軸力N (x)の厳密解は，える（既知の）値であることを意味している．

N (x ) = -P ( 0 £ x £ L ).

長柱の支持条件の数学モデルは，これらの境界条件式を組み合わせて，次のように理想化される．

式(17)を式(16)に代入すれば，中心軸圧縮力の作用を受け

ピン支点（ヒンジ支点）：

w = 0 on w,

M = 0 on M.

る均質かつ等断面長柱の線形座屈問題の支配方程式は，次(12)のように表される．

固定端：

EI

w = 0 on w,

q = 0 on .

(17)

(13)

d 4w (x )d2w (x )+P= 0 in W .
4
dx
dx 2

(18)

また，境界条件式を面外変位wで表すと，
自由端：

M = 0 on M, Q = 0 on Q.

ピン支点（ヒンジ支点）：(14)

w = 0 on w,

これらは，斉次の境界条件式である．

d2w= 0 on M.
dx 2

(19)

dw
= 0 on .
dx

(20)

-EI

固定端：
(4) 一般化応力と変位で表された構成方程式
長柱に対して，本稿のように右手直交直線座標系O-xyz

w = 0 on w,を設定すると，yz断面に設定した図心軸まわりの断面一次モーメントは零になる．このことから，一般化応力と変位で表された構成方程式に含まれる軸変位u (x)と曲げたわみ曲線を表す面外変位w (x)に関する項は非連成化され，微小

自由端：
-EI

ひずみ・微小回転の有限変位に基づく梁－柱理論の構成方(21)

程式は，次のように表される．
2üìïï d u (x ) 1 é d w (x ) ù ïï,ú ïN (x ) = EA ïí
+ êê
ú ý
ï
ï
d
x
2
d
x
ï
ï
ë
û ï
ï
î
þ

d2w
= 0 on M ,
dx 2

3

-EI

(15a)

dw
dw
-P
= 0 on Q.
3dxdx

本問題は，定数係数の四階線型の斉次常微分方程式の斉―28―

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

第 56 号(令和元年 11 月)

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

ここで，w0, wLは，それぞれ，長柱の左端点(x = 0)と右端求められた中心軸圧縮力を座屈荷重Pcrと呼び，これは長柱点(x = L)の面外変位wの境界条件に関するLagrangeの未定が安定平衡を与える中心軸圧縮力の上限を意味する．なお， 乗数であり，0, Lは，それぞれ，長柱の左端点(x = 0)と座屈荷重Pcrは固有値問題の固有値であり，これに対応する

右端点(x = L)のたわみ角の境界条件に関するLagrangeの固有関数w (x)が曲げ座屈変形の形状を表す．

未定乗数である．式(27)の右辺のバーを付した面外変位wとたわみ角は境界での既知の値であり，その下付き文字の0とLは，それぞれ，x = 0, x = Lでの条件であることを意

４．強形式と等価な変分原理の汎関数味している．斉次の基本境界条件を考えるならば，
P* = U + V(1) 許容な試行関数を仮定した変分原理の汎関数
均質かつ一様断面の長柱が中心軸圧縮力Pの作用により+lw 0 w

僅かに曲がったときの汎関数（全ポテンシャルエネルギ）+lq 0 q

x =0

x =0

+ lwL w

+ lqL q

x =L

.

(28)

x =L

は，次のように与えられる．

P = U +V .

(22)(3) 仮想ばね法により基本境界条件を緩和させた変分原理の汎関数

ここで，Uは長柱に蓄積されるひずみエネルギであり，Vは

仮想ばね法（penalty関数法）では，Lagrangeの未定乗数

中心軸圧縮力Pによる外力ポテンシャルである．ひずみエを次のように近似する22)．

ネルギUは，次のように表される．
L

l @ ag (v ) .

2

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

(23)

0

中心軸圧縮力Pの作用により長材が僅かに曲げ座屈変形を生じたときの右端部(x = L)の軸圧縮力の作用方向の圧縮
変位uLは，(29)

ここで，g (v) = 0は拘束条件であり，  は正の大きな値（penalty数）である．数学的には → ∞ のとき，g (v) → 0
であるが，数値計算ではに正の大きな値を設定したとき，
g (v) @ 0になることを期待する．本稿では，このpenalty数
を長柱の左端点(x = 0)と右端点(x = L)での非斉次の基本境

L界条件に対応する仮想のばね定数として考える17)．このと2

1 ó æ dw ö
DuL = - ô ççç ÷÷÷ dx .
÷2ôõ çè dx ø(24)

きの汎関数 P は，次のように表される．

0

P = U +Vaa.
+ w 0 (w - w 0 )2
+ wL (w - w L )2
x =0
x =L
2
2
b
b
+ q 0 (q - q0 )2+ qL (q - qL )2
x =0
x =L
2
2

よって，中心軸圧縮力Pによる外力ポテンシャルVは，
L

2P ó æ dw ö
V = -(-P )uL = - ô ççç ÷÷÷ dx .
÷2ôõ çè dx ø

(25)

0

(30)

式(23)と式(25)を式(22)に代入すれば，汎関数は次のよう

ここで，w0, wLは，それぞれ，長柱の左端点(x = 0)と右端に表される．

点(x = L)の面外変位wに抵抗する仮想の鉛直ばね定数であLL

2

り，0, Lは，それぞれ，長柱の左端点(x = 0)と右端点(x =2

EI ó æç d2w ö÷÷P ó æ dw ö÷ô çç
dx - ô çççP=÷÷÷÷ dx .
2÷2 ô2ôõ èç dx øõ èç dx ø÷(26)

次の基本境界条件を考えるならば，00L)のたわみ角に抵抗する仮想の回転ばね定数である．斉P = U +Vaa.
+ w0 w2+ wL w 2
=
x
x =L
0
2
2
b
b
+ q0 q2+ qL q 2
x =0
x =L
2
2

(2) Lagrangeの未定乗数法による付帯条件付き変分原理の汎関数
面外変位w (x)の近似関数に基本境界条件を満足しない
許容でない関数を仮定したとき，長柱の左端点(x = 0)と右

(31)

端点(x = L)での非斉次の基本境界条件を付帯条件とする汎関数*は，次のように表される．

P * = U +V
+lw 0 (w - w 0 )+lq 0 (q - q0 )

x =0

x =0

+ lwL (w - wL )+ lqL (q - qL )

x =L

.

(27)

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

x =L

式(31)の汎関数の変関数（長柱の面外変位）w (x)を次の―29―

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

第 56 号(令和元年 11 月)

試行関数W (x)で近似する．

由度で近似される．この近似された汎関数 P の停留条件：

w(x ) @ W (x )¶P
= 0 (i = 1,2, , nCP - 1, nCP ) ,
¶Ai

nCP

= å AN(x ) .
i i, p

(32)

(37)

i =1

= Ni , p A

より，次の離散化された座屈方程式を得る．

ここで，[(K0 + KP ) - PKG ]A = 0 .


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

A = {A1, A2 , , An -1, An }T ,
CP

(33)
ここで，

L

2

T

2

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

(34)

CP

(38)

K

であり，Tは転置を意味する．Ni, p (x)は次数pのB-spline基底(39)関数，iは基底関数の番号，Aiはi番目の基底の重み係数であは自由境界を有する長柱の剛性行列（すなわち，無拘束のり，nCP (= M + p + 1)は基底関数の数である．領域に1の分

長柱の剛性行列），割(partition of unity)を満たすB-spline基底を張るために必要

KP = aw 0 NiT, p N j , p

なノットの数Mは，
M = m + 2p .

+bq 0

(35)

x =0

dNiT, p dN j , pdxdx

付加ノットの数である．図-1を参照すれば，パラメータ領

x =L

,

dNiT, p dN j , p
dx

dx

(40)

x =L

は基本境界条件を緩和するための仮想ばねの剛性行列であり，K = K0 + KPは長柱の剛性行列である．また，

域[0, 1]の左側と右側に，それぞれ，p個の付加ノットが設

LT

ó dNi , p dN j , p
KG = P ôdx ,dx
õ 0 dx

定されていることが確認できよう．基底関数の数nCPと境界
を含む領域内に設定するノットの数mとの間には，次の関係が成り立つ．

+ bqL
x =0

ただし，mは領域内（境界を含む）に設定するノットの数であり，2pは領域の外側（境界を含まない）に設ける+ awL NiT, p N j , p(41)

は長柱の幾何剛性行列であり，0は零の列ベクトルである．
(36)

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

式(32)を式(31)の汎関数 P に代入すると， P は有限の自式(38)の自明でない解は，次式の座屈条件式（固有方程nCP = m + p - 1 .

Legendreの数値積分により求めた．

x, u

x, u
x, u

P

x, u

P

P

P

L
L

O

L

L

z, w
z, w

(a) P-P

z, w

z, w

O

O

O
(b) C-F(c) C-C図-3 数値実験に用いる長柱の支持条件―30―(d) C-P

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

第 56 号(令和元年 11 月)式）：

L

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

(42)

Pcr =

から固有値Pを求め，これを式(38)に代入して固有値Pに対応する固有ベクトルAを決定することで得られる．なお，

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

L

2

P ó æç dw ö÷ô çç ÷÷ dx
2ô
õ çè dx ÷ø

.

(46)

0

固有値と固有ベクトルのペアはnCP個だけ求められる．

ここで，分子は座屈による曲げ変形により長柱に貯えられるひずみエネルギ（内力仕事）であり，分母は長柱の曲げ
座屈変形により生じる軸圧縮変位と中心軸圧縮力とがな

６．数値実験および考察

す外力仕事を意味する．式(46)の右辺の表示は，線形化さ線形座屈問題を本手法により解析し，座屈荷重Pcrに与えるれた梁－柱理論に基づく長柱の線形座屈問題におけるRayleigh商である24)．この式から次のことがわかる．

spline次数p，境界を含む領域内に設定するノットの数（以

ここでは，中心軸圧縮力を受ける均質かつ等断面長柱の強形式を満足する面外変位w (x)の厳密解を用いれば，下，ノットの数と略す）m，仮想ばね定数の値および長柱式(46)により座屈荷重Pcrの厳密解を得ることができの支持条件の影響を詳細に調べる．

る24)．

本稿では長柱の座屈荷重Pcr ，仮想の鉛直ばね定数 w0,
wLおよび仮想の回転ばね定数0, Lを，次のように無次元化して表している．

なお，許容である試行関数を仮定したRitz法による近似解
PcrL2

P* =

EI

座屈荷重Pcrは，外力仕事に対するひずみエネルギの大きさとして評価される．

は，厳密解の上界を与える8), 10)ので，常に次の関係が成り
(43)

.

立つ．

Pcr ³ PcrExact .
3

a* =

b* =

aw 0L
EI

bq 0L

EI

(47)

3

=

=

awLL

,

EI
b qL L

EI(44)

ここで，不等式の左辺はRitz法による近似解であり，右辺は厳密解である．
本手法により求められる座屈荷重Pcrは，式(46)の右辺の

.

(45)

分母に対する分子の値の大きさにより決定される．結局のところ，式(38)を数値的に解いて得られる未定係数Aの値

数値実験に設定する支持条件は，図-3に示す工学的に重が座屈荷重Pcrの精度を決めていることになる．単純ノット

要な四つのケースとし，例えばC-Pのように表す．これはベクトルtからなるB-spline基底関数は領域の境界で有柱の下端(x = 0)で固定Cであり，上端(x = L)でピン支点（ヒ限の値を持つ基底関数がp個だけ存在することは，先に述

ンジ支点）Pであることを意味する．また，自由端はFで表べた通りである（図-1を参照）．例えば，p = 2次のB-spline

す．なお，(a)と(c)は対称の支持条件，(b)と(d)は非対称の境界条件であり，(a)と(b)は静定，(c)と(d)は不静定である．

基底関数をRitz法の近似解に仮定すると，長柱の下端境界(x = 0)を定義域とする基底は，一番目の基底N1,2 (x)と二番

通常，工学分野の諸問題の解決を目的とするならば，強目の基底N2,2 (x)である．正確には三番目の基底N3,2 (x)も境形式の近似解は有効数字三桁から四桁程度の精度があれ

界(x = 0)を定義域に持つが，N3,2 (0) = 0であるので今回は考ば十分である．しかし，渡辺23)が述べているように，Ritz法

える必要がない．故に，w (0) = 0の基本境界条件を満足すにより数値的に安定して高精度な近似解が計算できるよるためには，次式が成立しなければならない．

うになれば，Ritz解は厳密解やこれの級数形式の解に代わW (0) = A1 N 1,2 (0) + A2 N 2,2 (0) = 0 .



って有限要素の性能評価や各種離散化手法が強形式を正=1

しく解いているかの検証(verification)への活用が期待され

2

=1

(48)

2

る．このような視点から，本稿では無次元座屈荷重P*を有

本手法は仮想ばね法により基本境界条件を緩和させてい効桁数十桁程度で整理した．なお，数値計算には64 bit OS

るに過ぎないので，式(48)は近似的に満足されることになのpersonal computerを用い，Fortranの倍精度で計算した．

る．すなわち，式(31)の汎関数が停留するように求められた未定係数Aの成分A1とA2を用いて，式(48)が数値的に零に(1) Timoshenkoのエネルギ法に基づく本手法の近似解の精度評価
Timoshenokoのエネルギ法によれば，中心軸圧縮力を受ける均質かつ等断面長柱の座屈荷重Pcrは，次式によって求めることができる24)．

なることを期待する．近似解の精度を改善するために
spline次数pを高めると式(48)の項数は増加し，図-1からも分かるように基底関数Ni,p (x)の最大値は小さくなる．このことから，spline次数pを高次にすると数値計算上の誤差による近似解の精度低下が予想される．

―31―

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

第 56 号(令和元年 11 月)

表-1 中心軸圧縮力を受ける長柱の無次元座屈荷重Pに与えるspline次数p, ノット数のmと仮想ばね定数の影響: P-P
p



m

Mode 1 (S1)Mode 2 (A1)

Mode 3 (S2)Mode 4 (A2)

4

108

11

9.86960468

39.4784999

88.8290641

157.948063

21

9.86960441

39.4784188

88.8264727

157.914029

31

9.86960440

39.4784177

88.8264424

157.913700

41

9.86960440

39.4784176

88.8264401

157.913675

51

9.86960440

39.4784176

88.8264397

157.913672

1010

1012

5

8

10

10

10

12

10

61

9.86960440

39.4784176

88.8264397

157.913671

11

9.86960468

39.4784999

88.8290641

157.948063

21

9.86960441

39.4784188

88.8264727

157.914029

31

9.86960440

39.4784177

88.8264424

157.913700

41

9.86960440

39.4784176

88.8264401

157.913675

51

9.86960440

39.4784176

88.8264397

157.913672

61

9.86960440

39.4784176

88.8264397

157.913671

11

9.86960424

39.4784995

88.8290636

157.948063

21

9.86960433

39.4784187

88.8264726

157.914029

31

9.86960437

39.4784177

88.8264424

157.913700

41

9.86960437

39.4784176

88.8264401

157.913675

51

9.86960440

39.4784176

88.8264397

157.913672

61

9.86960439

39.4784176

88.8264396

157.913671

39.4784187

88.8265251

157.915905

11

9.86960440

21

9.86960440

39.4784176

88.8264398

157.913675

31

9.86960440

39.4784176

88.8264396

157.913671

41

9.86960440

39.4784176

88.8264396

157.913670

51

9.86960440

39.4784176

88.8264396

157.913670

61

9.86960440

39.4784176

88.8264396

157.913670

11

9.86960440

39.4784187

88.8265251

157.915905

21

9.86960440

39.4784176

88.8264398

157.913675

31

9.86960440

39.4784176

88.8264396

157.913671

41

9.86960440

39.4784176

88.8264396

157.913670

51

9.86960440

39.4784176

88.8264396

157.913670

61

9.86960440

39.4784176

88.8264396

157.913670

11

9.86960446

39.4784187

88.8265252

157.915905

21

9.86960439

39.4784176

88.8264398

157.913675

31

9.86960442

39.4784176

88.8264396

157.913671

41

9.86960440

39.4784176

88.8264396

157.913670

51

9.86960439

39.4784176

88.8264396

157.913670

61

9.86960440

39.4784176

88.8264396

157.913670

9.86960440

39.4784176

88.8264396

157.913670

Exact

25)

(2) 静定の長柱するのは最低次の座屈荷重であるが，二次以上の座屈荷重

まず，面外変位wに対応する無次元仮想鉛直ばね定数

*とこれに対応する座屈モードは，座屈荷重を高めるためのと無次元座屈荷重P，離散化条件の関係を整理する．表-1

重要な情報を提供してくれる．このことから，本稿では一は，支持条件がP-Pである中心軸圧縮力の作用による長柱の一次から四次までの無次元座屈荷重Pに与えるspline次

次から四次までの座屈荷重について検討した．なお，表中のS1とS2は，それぞれ，長柱の中央点(x = L/2)に関して対数p, ノットの数mと無次元仮想鉛直ばね定数 * の影響を称な一次と二次の座屈モードを，A1とA2は，それぞれ，同テーブル形式で整理したものである．分岐座屈荷重に相当点に関して逆対称な一次と二次の座屈モードを意味する．

―32―

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

第 56 号(令和元年 11 月)

表-2 中心軸圧縮力を受ける長柱の無次元座屈荷重Pに与えるspline次数p, ノットの数mと仮想ばね定数の影響: C-F
p



m

Mode 1Mode 2Mode 3Mode 44108

112.46740105

22.2066177

61.6856094

120.913496

21

2.46740105

22.2066096

61.6850340

120.902774

31

2.46740105

22.2066095

61.6850269

120.902662

41

2.46740105

22.2066095

61.6850264

120.902653

51

2.46740105

22.2066095

61.6850263

120.902652

1010

1012

5

8

10

10

10

12

10

61

2.46740105

22.2066095

61.6850263

120.902652

11

2.46740110

22.2066181

61.6856105

120.913499

212.46740110

22.2066100

61.6850352

120.902776

31

2.46740110

22.2066099

61.6850282

120.902664

41

2.46740110

22.2066099

61.6850276

120.902656

51

2.46740110

22.2066099

61.6850275

120.902654

61

2.46740110

22.2066099

61.6850275

120.902654

11

2.46740126

22.2066196

61.6856146

120.913507

212.46740102

22.2066093

61.6850333

120.902772

31

2.46740111

22.2066100

61.6850283

120.902664

41

2.46740106

22.2066095

61.6850265

120.902654

51

2.46740111

22.2066100

61.6850278

120.902655

61

2.46740109

22.2066099

61.6850274

120.902654

61.6850370

120.903095

11

2.46740105

22.2066095

21

2.46740105

22.2066095

61.6850263

120.902653

31

2.46740105

22.2066095

61.6850263

120.902652

41

2.46740105

22.2066095

61.6850263

120.902651

51

2.46740105

22.2066095

61.6850263

120.902651

61

2.46740105

22.2066095

61.6850263

120.902651

112.46740110

22.2066100

61.6850383

120.903098

21

2.46740110

22.2066099

61.6850275

120.902655

31

2.46740110

22.2066099

61.6850275

120.902654

41

2.46740110

22.2066099

61.6850275

120.902654

51

2.46740110

22.2066099

61.6850275

120.902654

61

2.46740110

22.2066099

61.6850275

120.902654

11

2.46740135

22.2066121

61.6850440

120.903108

212.46740114

22.2066103

61.6850287

120.902657

31

2.46740105

22.2066095

61.6850262

120.902651

41

2.46740109

22.2066098

61.6850273

120.902654

51

2.46740109

22.2066098

61.6850272

120.902653

61

2.46740111

22.2066100

61.6850277

120.902654

Exact25)

2.46740110

22.2066099

61.6850275

120.902654

無次元仮想鉛直ばね定数*は，108, 109, 1010, 1011, 1012の五
8

10

12

関数まで連続になるので，せん断力Q (x)は一次の区分的多ケースについて検討した．本稿では，10 , 10 , 10 の結果

項式で近似される．当然のことではあるが，spline次数p =を示している．梁－柱理論に基づく線部材の二節点有限要5に設定すると，面外変位wは四階導関数まで連続であるの素に用いられるHermite型の形状関数は，隣接する要素間でで，せん断力Q (x)の近似式は二次の区分的多項式になる．

面外変位wとたわみ角が連続になるような三次の区分的ノットの数mは，11から61まで変化させ，有効桁数九桁で多項式である．このことを鑑み，本稿ではspline次数pを4と表した無次元座屈荷重Pの値の変化を調べた．表のboldの5に設定した．Spline次数p = 4のとき，面外変位wは三階導

数値は近似解が収束し，かつ厳密解と一致したことを意味―33―大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

表-3 中心軸圧縮力を受ける長柱の無次元座屈荷重Pに与えるspline次数p, ノットの数mと仮想ばね定数の影響: C-C
p



m

Mode 1 (S1)Mode 2 (A1)

Mode 3 (S2)Mode 4 (A2)

4

108

11

39.4785116

80.7649907

157.951191

238.965477

2139.4784173

80.7629243

157.914049

238.720260

31

39.4784161

80.7629003

157.913695

238.718188

41

39.4784160

80.7628985

157.913669

238.718045

51

39.4784160

80.7628982

157.913665

238.718023

1010

12

10

5

8

10

10

10

12

10

61

39.4784160

80.7628981

157.913665

238.718018

1139.4785131

80.7650066

157.951197

238.965525

21

39.4784189

80.7629403

157.914055

238.720307

3139.4784177

80.7629163

157.913701

238.718235

41

39.4784176

80.7629144

157.913676

238.718092

51

39.4784176

80.7629142

157.913672

238.718071

61

39.4784176

80.7629141

157.913671

238.718066

11

39.4785182

80.7650200

157.951218

238.965574

21

39.4784168

80.7629363

157.914047

238.720295

31

39.4784186

80.7629182

157.913738

238.718241

41

39.4784169

80.7629132

157.913673

238.718089

51

39.4784179

80.7629150

157.913673

238.718073

61

39.4784174

80.7629138

157.913670

238.718065

11

39.4784169

80.7629416

157.915538

238.740313

2139.4784160

80.7628982

157.913668

238.718053

31

39.4784160

80.7628981

157.913664

238.718017

41

39.4784160

80.7628981

157.913664

238.718016

51

39.4784160

80.7628981

157.913664

238.718016

61

39.4784160

80.7628981

157.913664

238.718016

1139.4784186

80.7629579

157.915545

238.740361

2139.4784176

80.7629142

157.913674

238.718100

3139.4784176

80.7629140

157.913670

238.718064

41

39.4784176

80.7629141

157.913670

238.718063

51

39.4784176

80.7629141

157.913670

238.718063

61

39.4784176

80.7629141

157.913670

238.718063

11

39.4784358

80.7629941

157.915616

238.740470

21

39.4784191

80.7629175

157.913681

238.718111

31

39.4784157

80.7629104

157.913663

238.718053

41

39.4784174

80.7629138

157.913670

238.718063

51

39.4784172

80.7629134

157.913669

238.718061

61

39.4784180

80.7629151

157.913672

238.718066

39.4784176

80.7629142

157.913670

238.718064

Exact

25)

する．なお，厳密解25)は各支持条件に対応する座屈条件式

鉛直ばね定数* = 108, 1010である．このとき，座屈荷重のを数値的に解き，有効桁数十五桁で収束した座屈荷重の値

近似解は上界であり，ノットの数mの増大にともなって上を求めた．厳密解の値もboldで示してある．この支持条件

から厳密解に近づいていく．同じspline次数の無次元仮想の曲げ座屈モードの厳密解は，正弦関数による弾性変形項

ばね定数* = 1012のときの一次の座屈荷重は一定の値に収のみで表される．

束せずに振動しており，m = 21, 51の一次の座屈荷重は下界

表-1より，一次から四次までの座屈荷重の最も良好な収値を与えている．これは，p = 4のときも同様である．これ

束状態を示しているのは，spline次数p = 5かつ無次元仮想より，無次元仮想鉛直ばね定数* = 1012に設定するとsplime―34―大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

表-4 中心軸圧縮力を受ける長柱の無次元座屈荷重Pに与えるspline次数p, ノットの数mと仮想ばね定数の影響: C-P
p



m

Mode 1Mode 2Mode 3Mode 44108

11

20.1907334

59.6800317

118.910020

197.959034

21

20.1907274

59.6795193

118.899978

197.858764

31

20.1907274

59.6795130

118.899872

197.857875

1010

1012

5

8

10

10

10

12

10

41

20.1907273

59.6795125

118.899864

197.857812

51

20.1907273

59.6795124

118.899862

197.857803

61

20.1907273

59.6795124

118.899862

197.857800

1120.1907346

59.6800351

118.910027

197.959045

21

20.1907286

59.6795229

118.899985

197.858776

31

20.1907286

59.6795165

118.899879

197.857887

41

20.1907285

59.6795160

118.899871

197.857824

51

20.1907285

59.6795159

118.899870

197.857814

61

20.1907285

59.6795159

118.899869

197.857812

11

20.1907355

59.6800390

118.910035

197.959060

21

20.1907282

59.6795211

118.899981

197.858770

31

20.1907285

59.6795167

118.899879

197.857888

41

20.1907282

59.6795150

118.899869

197.857821

51

20.1907287

59.6795162

118.899870

197.857815

61

20.1907285

59.6795158

118.899869

197.857812

11

20.1907274

59.6795214

118.900268

197.865558

21

20.1907273

59.6795124

118.899863

197.857813

31

20.1907273

59.6795124

118.899862

197.857800

41

20.1907273

59.6795124

118.899862

197.857799

51

20.1907273

59.6795124

118.899862

197.857799

61

20.1907273

59.6795124

118.899862

197.857799

1120.1907286

59.6795250

118.900275

197.865570

21

20.1907286

59.6795160

118.899870

197.857825

31

20.1907285

59.6795159

118.899869

197.857811

41

20.1907285

59.6795159

118.899869

197.857811

51

20.1907285

59.6795159

118.899869

197.857811

61

20.1907285

59.6795159

118.899869

197.857811

11

20.1907306

59.6795306

118.900285

197.865584

21

20.1907290

59.6795171

118.899873

197.857829

31

20.1907281

59.6795147

118.899867

197.857808

41

20.1907285

59.6795158

118.899869

197.857811

51

20.1907285

59.6795157

118.899869

197.857810

61

20.1907286

59.6795162

118.899870

197.857812

20.1907286

59.6795159

118.899869

197.857811

Exact

25)

次数pにかかわらず，若干ではあるが数値解に影響を与え条件がC-Fである長柱の一次から四次までの無次元座屈荷ていると考えられる．なお，spline次数p = 4のときの四次の重Pに与えるspline次数p, ノットの数mと無次元仮想ばね

座屈荷重は無次元仮想鉛直ばね定数 の値にかかわらず，

定数*の影響を表-2に示す．ここでは* = *と設定し，そノットの数mを増やしても一定値に収束しない．以上の考れを*で表すこととする．表-1の結果から判断すると，無察は，あくまで有効桁数を九桁としたときのものである．

次元仮想鉛直ばね定数 * は108 と1010 の二ケースについて

次に，たわみ角に対応する無次元仮想回転ばね定数*

検討すればよいであろうが，参考のために1012の結果も示*

と離散化条件，無次元座屈荷重P の関係を整理する．支持した．なお，表-1と同様に，boldの数値は近似解が収束し，―35―

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

第 56 号(令和元年 11 月)

かつ厳密解と一致したことを意味している．ただし，近似

表-4は，表-3の上端のたわみ角の拘束を解除したモデル解は収束したが，その収束値が厳密解と一致しないものにに相当する．この支持条件の曲げ座屈モードの厳密解も余は，boldの数値に下線を付した．この支持条件の曲げ座屈弦関数と正弦関数による二つの弾性変形，剛体回転と剛体モードの厳密解は，余弦関数による弾性変形と剛体変位の変位の四つの項で表される．この支持条件についても
spline次数p = 5かつ無次元仮想ばね定数* = 1010の組み合わ二つの項により表現される．
表-2より，無次元仮想ばね定数* = 1010のとき，spline次せが最も良好な結果を与えている．表-3とは異なり，二次数pに係わらず，一次から四次までの座屈荷重は上界であから四次の座屈荷重の収束値は上界であるが，一次のそれり，ノットの数mの増大にともなって上から厳密解に近づは下界である．ただし，これも有効桁数八桁で整理すれば，

く良好な収束状態を示している．その収束値は，有効桁数

厳密解に一致する精度の高い上界値になる．

九桁で厳密解と一致する．他方，無次元仮想ばね定数* =(2)と(3)より，無次元仮想ばね定数* = 1010に設定し，離108のときは，仮想ばねによる長柱の剛性行列の評価が正散化条件をspline次数p = 5かつノットの数m = 41に設定ししくなく，ノット数m = 11の無次元座屈荷重は上界であるたときの長柱の一次から四次までの無次元座屈荷重P*は，が，座屈荷重は厳密解より小さな下界の値に収束している． 支持条件にかかわらず有効桁数九桁での収束値である．そこれは，表-1と異なる点である．なお，無次元仮想鉛直ばの収束値は，有効桁数八桁から九桁の範囲で厳密解に一致

ね定数 を10 と10 に設定し，それぞれに対して無次元仮するという高い近似精度が示された．

*

8

10

想回転ばね定数 * を108 から1012 まで変化させてみたが結論は変わらなかった．

７．結言(3) 不静定の長柱(2)より，静定の長柱の座屈荷重の近似解に影響を与えな

本稿では，線形化された梁－柱理論に従う中心軸圧縮力

い無次元仮想鉛直ばね定数*の最大値は1010であり，固定を受ける均質かつ真っすぐな長柱の線形座屈問題を例題

端については無次元仮想鉛直ばね定数*，無次元仮想回転とし，この強形式の厳密解を正解として仮想ばね法により基本境界条件を緩和させたB-spline Ritz法の座屈荷重の精

ばね定数*共に1010であることがわかった．ここでは，この結果に基づき不静定な長柱について(2)と同様の検討を行

度，spline次数，解析領域内に配置するノットの数，仮想ばう．表-3と表-4は，それぞれ，C-C, C-Pの支持条件を有するね定数の値および支持条件の関係を詳細に調べた．本調査

長柱の一次から四次までの無次元座屈荷重P に与えるから，次の知見が得られた．
[1] 支持条件とspline次数pにかかわらず，数値的に安定しspline次数p, ノットの数mと無次元仮想ばね定数 *の影響て近似解に影響を与えない範囲の仮想ばね定数の最を整理したものである．

大値は，無次元仮想鉛直ばね定数*，無次元仮想回転

表-3は表-2と似たような傾向を示しており，最も良好な

ばね定数*共に1010であった．

結果を与えるのは，spline次数p = 5かつ無次元仮想ばね定数* = 1010の組み合わせであろう．この支持条件の一次と[2]三次の対称な座屈モードの正解は表-2と同様の二項によSpline次数p = 5かつ無次元仮想ばね定数* = 1010に設定したときの長柱の一次から四次までの無次元座屈

り表現され，無次元座屈荷重P は有効数字九桁で収束し，

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

その値は厳密解に一致する．これに対して，二次と四次の大にともなって上から厳密解に近づくような収束状

逆対称の座屈モードは余弦関数と正弦関数による二つの態を示し，有効桁数九桁での収束値が得られる．ただ弾性変形項，剛体回転項と剛体変位項の四つの混ぜ合わせし，この収束値は厳密解の下界の値を与えることがにより表現され，対称座屈モードに比べてやや複雑である．

ある．
[2]で得られた収束値は，有効桁数八桁から九桁の範

このときの座屈荷重も有効数字九桁で収束値を得ている[3]

が，その値は厳密解よりも僅かに小さい．表-1，表-2と表囲で本問題の座屈荷重の厳密解に一致する．

-3の一次と三次の座屈荷重の収束状態と精度から判断す

なお，固有関数とその導関数の精度についても検討していると，二次と四次の座屈荷重が下界となる原因はその座屈るが，これについては別の機会に報告する．今後は，基本

モードの形状が影響していると考えられる．この考察も有境界条件を正確に満足する許容な試行関数を仮定したB-

効桁数を九桁としたときのものであって，有効桁数を八桁にすれば，spline次数p = 5かつ無次元仮想ばね定数* = 1010

spline Ritz法により同様の調査を実施し，本稿の結果と比較する予定である．
最後に，本稿で示した手法に新規性はないが，Ritz解をに設定した本手法の収束値は厳密解に一致する．なお，表-1のP-P, 表-2のC-Fと異なり，C-Cのspline次数p = 4かつ無厳密解や級数解の代替として活用する 23) という従来あま

次元仮想ばね定数* = 1010の二次から四次までの座屈荷重り検討されることのなかった視点から調査結果を整理しは有効桁数九桁での収束値が得られておらず，不静定の長た．本稿は同研究分野の専門家による査読を経ずに公開さ

柱は近似解の収束性が悪くなっていることも確認できる．

れるが，本稿で示した表の数値は，読者の判断で種々の離―36―大分工業高等専門学校紀要

第 56 号(令和元年 11 月)

12)

散化手法の検証などに使っていただければ幸いである．

H. Nagino, T. Mikami, T. Mizusawa: Three-dimensional
free vibration analysis of isotropic rectangular plates using

謝辞：本研究の一部は，LIXIL住生活財団調査研究助成

the B-spline Ritz method, Journal of Sound and Vibration,（2013年度）および前田記念工学振興財団研究助成（平成Vol.317 (1-2), pp.329-353, 2008.

27年度）を受けて行われました．ここに記して謝意を示し

13)

る中実円筒体の3次元自由振動解析，構造工学論文集,
Vol.54A, pp.90-101, 2008.

ます．
14)

参考文献
1)

2)

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

Using the Lagrange Multiplier Method, Journal of Soundand Vibration, Vol.173 (4), pp.457-470, 1994.
15)

16)

学研究会シンポジウム論文集 (USB), R4-1 (6 pages),
2017.
山本寧音，名木野晴暢，足立忠晴，樋口理宏：軸方向

17)
18)

6)

山本寧音，名木野晴暢，足立忠晴，岩崎英治：軸方向

plates by using B-spline functions, Journal of Sound andVibration, Vol.62 (2), pp.301-308, 1979.
20)

A.W. Leissa: The historical bases of the Rayleigh and Ritz

2005.

and mesh refinement, Computer Methods in Applied
Mechanics and Engineering, 194 (39–41), pp.4135-4195,
21)

Wiley & Sons, 2009.

菊地文雄：有限要素法概説 [新訂版] －理工学におけ22)

pp.185-197, 2000.
23)

図書出版, pp.76-77, 1979.

11)

T.J.R. Hughes: The Finite Element Method: Linear Staticand Dynamic Finite Element Analysis, Dover Publications,

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

10)

J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs: IsogeometricAnalysis: Toward Integration of CAD and FEA, Johnる基礎と応用－ (FEM + BEM = 3), サイエンス社,
pp.25-38, 2007.

9)

T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs: Isogeometricanalysis: CAD, finite elements, NURBS, exact geometry

pp.961-978, 2005.

8)

T. Mizusawa, T. Kajita, M. Naruoka: Vibration of skew

傾斜機能材料からなる線状部材の線形座屈解析，土木 学 会 第 73 回 年 次 学 術 講 演 会 講 演 概 要 集 (DVDROM), I-538, pp.1075-1076, 2018.
methods, Journal of Sound and Vibration, Vol.287 (4-5),
7)

R. Kao: Application of hill functions to two-dimensional
Structures, Vol.11 (1), pp.21-31, 1975.

19)

件の影響，土木学会西部支部 研究発表会講演概要集(CD-ROM), I-23, pp.45-46, 2018.
5)

用性, 応用力学論文集, Vol.12, pp.143-154, 2009.
R. Kao: Approximate solutions by utilizing hill functions,plate problems, International Journal of Solids and山本寧音，名木野晴暢，足立忠晴：軸方向傾斜機能材料からなる線状部材の線形座屈特性に与える支持条

名木野晴暢，水澤富作，三上隆：点支持された周辺自由Mindlin平板の自由振動解析へのB-spline Ritz法の適Computers & Structures, Vol.3 (2), pp.397-412, 1973.

に傾斜機能を有する線構造部材の線形座屈固有値問
題の解析手法，土木学会西部支部 研究発表会講演概要集 (CD-ROM), I-6, pp.11-12, 2017.

内山武司, 上田正生：点支持を有するMindlin平板の自由振動解析, 日本建築学会構造系論文集, 63巻, 512号,
pp.83-89, 1998.

山本寧音，名木野晴暢，足立忠晴，樋口理宏：軸方向る縦弾性係数の変化の影響，第5回九州橋梁・構造工

4)S. Kitipornchai, Y. Xiang, K.M. Liew: Vibration Analysis

土木学会構造工学委員会計算力学とその応用に関す

傾斜機能材料からなる柱部材の線形座屈特性に与え3)

名木野晴暢，三上隆，水澤富作：B-spline Ritz法によ渡辺力：p-Ritz法における試行関数の改良に関する一考察，構造工学論文集A, 61A巻, pp.1-10, 2015.

林毅，村外志夫：変分法 (応用数学講座 第13巻), コ

24)

ロナ社, pp.1-89, 1958.
C.L.ディム，I.H.シャームス (共著)，砂川恵 (監訳)：9), 技報堂出版, pp.37-98, 1982.
25) C.M. Wang, C.Y. Wang, J.N. Reddy: Exact solutions for

材料力学と変分法 (理工学海外名著シリーズ26), ブ

buckling of structural members, CRC Press, pp.9-75

レイン図書出版, pp.94-151, 1977.

(chapter 2), 2004.

山崎徳也，彦坂煕 (共著)：構造解析の基礎, 共立出版,
pp.165-187, 1978.

―37―

福本秀士：構造物の座屈・安定解析 (新体系土木工学(2019.9.30 受付)