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

第 47 号(平成 22 年 11 月)

外部で液体に接する中実円筒梁の固有円振動数について－ 液体の高さと梁の長さが異なる場合 －金城 和久1・名木野 晴暢2・水澤 富作3・三上 隆41琉球大学大学院，2本校都市システム工学科，3大同大学，4北海道大学大学院

橋梁の橋脚などに多く見られる外部で液体に接する円筒構造物の構造設計には，基本的な情報として，液体の影響を考慮した構造物の固有円振動数を正確に把握することが不可欠である．現実的には，液体の高さと構造物の長さが異なる問題を取り扱うことになる．この問題の厳密解は，液体接触領域と非接触領域の二つの基礎方程式を連続条件と境界条件の下で解けば求められるが，解式が複雑になり，実用的とは言い難い．本報告では，べき級数を試行関数に採用したRitz法 (p-Ritz法) を用いて，外部で任意の高さの液体に接する中実円筒梁の自由振動問題を定式化し，実用的な構造解析プログラムを開発した．
また，本解析法による数値解の収束状態は良好であり，高い解析精度が得られることを確認した．

キーワード : 流体－構造相関自由振動問題，中実円筒梁，固有円振動数，p-Ritz法

１．まえがきr, w



外部で液体に接する円筒構造物は，橋梁の橋脚や海洋構
造物などに用いられる．これらの構造物は，地震動や波浪Roなどの複雑な動的外力を受けるため，その構造設計では，基礎的な情報として，液体の影響を考慮した構造物の自由振動特性 (固有円振動数と固有振動モード) や構造物と液x体の動的相互作用で生じる動水圧を出来るだけ正確に把握することが重要な課題となる．
桜井1) は，水中にある梁の自由振動問題の解析理論を構L築し，その解析手法の例として，数学的に閉じた解 (厳密H解) を求める手順とエネルギー原理に基づくRayleigh法による解析手順について述べており，Rayleigh法を用いて，近似的に液体の影響を考慮した梁の固有円振動数を求め

oている．桜井1) の研究では，液体の高さと梁の長さが等し

図-1 解析モデルと円筒座標系

い場合の問題のみの取り扱いであるが，現実的には，液体の高さと梁の長さが異なる場合の問題を解明する必要がる検討を行った．本報告の結果は，有限要素法や種々の数あるであろう．しかし，この問題は，梁と液体が接触する
領域 (接触領域) と梁と液体が接触しない領域 (非接触領値解析法の数値解を評価するための参照解として活用できるように，全て表形式で示した．

域) で表される二つの基礎方程式を連続条件と境界条件の下で解かなければならないので，厳密解を導出する解式

２．液体の基礎方程式，境界条件および連成条件

はかなり複雑になり，実用的な解法とは言い難い．
本報告では，べき級数を試行関数に採用したRitz法(p-Ritz法) を用いて，任意の液体の高さを有する外部で液図-1には，解析モデルおよび円筒座標系 (r, , x) が示し体に接する中実円筒梁の自由振動問題を定式化し，実用的てある．ここで，L は円筒梁の長さ，Ro は外径であり，な構造解析プログラムを開発した．また，本解析法により

H は液体の高さである．下端 (x = 0) で固定，上端 (x = L)

求められる固有円振動数の収束状態やその妥当性に関すで自由な均質円筒梁の運動は調和振動するものとし，微小―30―

第 47 号大分工業高等専門学校紀要(平成 22 年 11 月)

変形かつ線形弾性であるBernoulli-Euler梁理論2) に従うものとする．また，円筒梁の面外変位成分を w (x, t) で表す．

ここで，次式で表される無次元座標系および無次元変位を導入する．

液体の運動は，非粘性，非回転かつ非圧縮性とすれば，
式(1) で表される速度ポテンシャル (r, , x, t) を仮定で

2

iw t

2

(1)

用いて，式(7) を無次元化すると，次式のようになる．

æ l r ö÷
ç
2K 1 çç j ÷÷÷æ l Lx ö÷¥çè H ÷ø
iw t
ç
F = iw e åcos q cos çç j ÷÷÷æöççè H ÷ø
l
R
j =0
÷
ç
lj K 1¢ çç j o ÷÷÷çè H ÷ø
æ l L x ö÷
H /L
ç
2
´ò
LW
(x)cos çç j ÷÷÷ dx
0
ççè H ÷ø2

¶ F 1 ¶F
1 ¶F ¶F
(2)
+
+
+
=0
¶r 2 r ¶r r 2 ¶q 2 ¶x 2
ここで， (r, , x) はポテンシャル，i2 = – 1 は虚数単位であり， は液体－円筒梁系の固有円振動数である．また，
動水圧 p (r, , x, t) は，次式で表される線形化された
Bernoulli の定理で与えられる．
p(r , q, x , t ) = -r f

¶F
¶t

(8)

W (x ) = LW (x )ただし，W ( ) は無次元振幅変位である．よって，式(8) をき，基礎方程式は式(2) で与えられる．
F(r , q, x , t ) = iw f(r , q, x )ex = x /L ,

(3)

(9)

ここで，H / L ≤ 1 は液高比であり，液体の高さを表す無次元量である．

ただし，f は液体の密度である．
解析モデルに対応する境界条件は，次式で与えられる．
(4)

３．p-Ritz法による構造－流体相関自由振動問題の定式化

ここで，式(4)は，それぞれ，液体の高さは十分に高く，重ここでは，p-Ritz法を用いて，構造－流体相関自由振動

p

x =H

= 0,

¶F
=0,
¶x x = 0

lim F = 0r ¥

力の影響を無視した自由表面条件式，x 軸方向の速度が液問題を定式化する．
液体－円筒梁系の全ポテンシャルエネルギー は，次底において零であるという液底での不透過条件式，振動数が高く，自由表面上の波が円筒梁の近傍から外側に伝播し式で与えられる．

ないという仮定に基づく液体外部接触問題に対する無限P = U max - (Tmax + Fmax )遠方放射条件式，に対応している．

(10)

円筒梁の r 方向の振動速度と外径 Ro における液体の速ここで，Umax は円筒梁の最大ひずみエネルギー，Tmax は度が等しいとすれば，円筒梁と液体との連成条件は，次の円筒梁の最大運動エネルギー，Fmax は液体の最大運動エ

ように与えられる．

ネルギーであり，これらは，次式で与えられる．

¶F
¶w
=
cos q
¶r r =R
¶t
r =R
o

2

(5)

o

ここで，
w(x , t ) = W (x )eU max =

iw t

(6)であり，W (x) は円筒梁の振幅変位である．
従って，基礎方程式の一般解を求め，境界条件と連成条

ì d2W (x )ï
ü
rs Aw 2L3 1 2
EI 1 ï
ï
ï
,
d
x
T
=
í
ý
ò 0 W (x)dx ,maxï
2
2L ò 0 ïïî dx ï
þæ l R ö÷çç j o ÷
K
2
2
÷÷ ì
çç
1
4
2 ¥
ü
æ l Lx ö÷ ïH /Lr L Ro pw
è H ÷ø ïïçç j ÷ ïïï
Fmax = f
W
x
x
(
)cos
d
÷ ý
å æl R ö ïíò 0çç
2
çè H ÷÷ø ïï
j =0
÷ ïï
ç
ï
þ
lj K1¢ çç j o ÷÷÷ îçè H ÷ø

件を満足するように解けば，液体と円筒梁の連成を考慮し(11)た速度ポテンシャル の閉じた形式の解析解が，次式のただし，EI, s, A は，それぞれ，梁の曲げ剛性，密度およように求められる．

び断面積であり，E は縦弾性係数，I は断面二次モーメン

æ l r ö÷
ç
2K1 çç j ÷÷÷æ l x ö÷ H
æ l x ö÷¥çè H ÷ø
iw tçç
cos q cos çç j ÷÷÷ ò W (x )cos çç j ÷÷÷ dx
F = iw e å0æöççç H ø÷
H
÷
l
R
ç
j =0
÷
è
ø
è
ç
lj K1¢ çç j o ÷÷÷çè H ÷ø

トである．
さて，Ritz法で仮定する変位関数は，少なくとも幾何学
的境界条件を満足した許容関数でなければならない．よって，幾何学的境界条件を満足する無次元振幅変位W ( ) は，
次のようなべき級数 (単純な多項式) で仮定する．

(7)¥

W (x ) = å as x 2 (2x - 1)sただし，j = {(2j + 1) } / 2 は余弦関数の零根，Kn (jr / H)s =0

は n 次の第二種変形Bessel関数であり，K’n (jRo / H) は
r = Ro におけるKn (jr / H) の r に関する一階微分である．

ここで，as は未定係数である．

―31―

(12)

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

表-1

第 47 号(平成 22 年 11 月)

外部で液体に接する短い中実円筒梁の振動数パラメータ の収束性に与えるべき級数の打ち切り項数S と級数の打ち切り項数J の影響と精度比較：L / Ro = 2,

H/L=1


J

S

1st

2nd

3rd

4th

5th

11

5
6
7
8
9
10
11
12
13
14
5
6
7
8
9
10
11
12
13
14
5
6
7
8
9
10
11
12
13
14

3.29603.29603.29603.29603.29603.29603.29603.29603.29603.29603.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.29483.294720.47120.47020.47020.47020.47020.47020.47020.47020.47020.47020.462
20.462
20.462
20.462
20.462
20.462
20.462
20.462
20.462
20.462
20.462
20.462
20.461
20.461
20.461
20.461
20.461
20.461
20.461
20.461
20.46159.92958.54058.54058.52458.52458.52358.52358.52358.52358.52359.90058.51558.51558.49858.49858.49858.49858.49858.49858.49859.89958.51458.51458.49758.49758.49758.49758.49758.49758.49758.497123.48123.21116.39116.39116.19116.19116.19116.19116.19116.19123.41
123.14116.34116.34116.14116.14116.14116.14116.14116.14123.40123.14116.34116.34116.14116.14116.14116.14116.14116.14116.14551.38216.47216.05194.80194.80193.65
193.65
193.63193.63193.63550.36216.29215.89194.71
194.71
193.57193.57193.54193.54193.54550.30216.29215.89194.71
194.71
193.57193.56193.54193.54193.54193.54

51

101

Exact

の数として，(2n – 1) 次の多項式まで厳密な積分値を求め

式(12) を式(11) に代入すれば，EI
U max ={a }T [K]{a } ,2LFmax =

r f L4Ro pw 2
2

Tmax =

rs Aw 2L3
2

ることができる．しかし，付加質量マトリックス[MF] は，
T

被積分関数がべき級数と級数項j に依存する余弦関数の{a } [M]{a } ,積の形になっているため，積分の厳密値を得るために必要な積分点の数n の算定式が存在しない．そこで，本報告でT{a } [MF ]{a }(13)点のGauss-Legendre の数値積分を用いて，積分値を求めることとした．当然の事ながら，この数値積分における積分を得る．ただし，{a }T = {a1 a2 a 3  aS }は，積分領域を幾つかの小区間に分割し，その区間毎に24(14)であり，S はべき級数の打ち切り項数，[K] は梁の剛性マ値の精度は，積分領域を分割する数に大きく依存する．
よって，式(13) を式(10) に代入し，未定係数ベクトル{a}T で極値化すれば，

トリックス，[M] は梁の質量マトリックスであり，[MF] は¶

液体の付加質量マトリックスである．全体基底関数を試行関数に採用した p-Ritz法より定式化される剛性マトリッ¶{a }T= {0}クス，質量マトリックスおよび付加質量マトリックスは，であり，次式のような線形代数方程式を得る．

対称の密行列になる．ここで，剛性マトリックス[K] と質éæ r ö÷æ ö
ïì
ïüïùú
ê
çç f ÷ç L ÷÷2ïïç÷
ê[K] - W í[M] + ç ÷ç ÷÷[MF ]ïýú {a } = {0}çç r ÷çïïïïúêè s ø÷çè Ro ÷ø
îï
þïûú
ëê

量マトリックス[M] は，被積分関数がべき級数のみで表されるため，Gauss-Legendre の数値積分により n を積分点

―32―

(15)

(16)

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

表-2

第 47 号(平成 22 年 11 月)

外部で液体に接する長い中実円筒梁の振動数パラメータ の収束性に与えるべき級数の打ち切り項数S と級数の打ち切り項数J の影響と精度比較：L / Ro = 40,

H/L=1


J

S

1st

2nd

3rd

4th

5th

51

5
6
7
8
9
10
11
12
13
14
5
6
7
8
9
10
11
12
13
14
5
6
7
8
9
10
11
12
13
14

2.98272.98272.98272.98272.98272.98272.98272.98272.98272.98272.98212.98212.98212.98212.98212.98212.98212.98212.98212.98212.98192.98192.98192.98192.98192.98192.98192.98192.98192.98192.981918.69518.69518.69418.69418.69418.69418.69418.69418.69418.69418.691
18.69018.69018.69018.69018.69018.69018.69018.69018.69018.69018.69018.68918.68918.68918.68918.68918.68918.68918.68918.68953.85152.51552.51452.49852.49852.49852.49852.49852.49852.49853.83852.50352.50252.48652.48652.48652.48652.48652.48652.48653.83552.50152.50052.48452.48452.48352.48352.48352.48352.48352.483110.12109.89103.51103.51103.32103.32103.32103.32103.32103.32110.08109.86103.49103.49103.30103.30103.30103.30103.30103.30110.08109.85103.49103.49103.29103.29103.29103.29103.29103.29103.29497.78
192.73
192.31172.75172.75
171.66
171.66
171.64171.64171.64496.96
192.65192.24172.71172.71
171.62
171.62
171.60171.60171.60496.79192.64
192.23172.70172.70171.61
171.61
171.59171.59171.59171.59

101

151

Exact

ここで，

本報告では，土木構造物の耐震設計を想定し，基本振動W2 = w 2L4

rs A
EI

(17)

から五次振動までの固有円振動数を対象としている．また，円筒梁の形状に関するパラメータを長さ－外径比L / Roと定義する．例えば，L / Ro = 2 は太く短い中実円筒梁，L / Roである．従って，式(18) の振動数方程式を解けば，固有値= 40 は細く長い中実円筒梁を意味する．

および固有ベクトルが求められる．

éìüù
æ r ÷öæ öïïïúêçç f ÷ç L ÷÷2ïïç÷
det ê[K] - W í[M] + ç ÷ç ÷÷[MF ]ïýú = 0
çç r ÷ç
ï
ïú
ê
R
÷
ç
è s ø÷è o ø
ïïûú
îïï
þ
ëê

数値計算では，コンクリート構造物を想定し，円筒体の密度s = 2.30 × 103 kg / m3，液体の密度f = 1.00 × 103 kg /(18)

m3 (密度比f / s は0.435) を用いる．なお，固有円振動数は，次式で表される振動数パラメータ で表す．

ここで，固有値計算には，Householder-QR法を用いており，各固有値に対応する固有ベクトルの計算には逆反復法を

W = wL2用いている．

rs A
EI

(14)

ここで，先にも述べたように，本報告で採用している数値
積分における積分値の精度は，積分領域を分割する数に大４．数値計算きく依存する．そこで，式(11) で表される液体の最大運動エネルギーFmax の被積分関数を対象として，種々の検討

ここでは，本解析法による数値解の収束状態およびその妥当性について検討する．

を行った結果，積分領域の分割数を15 に設定すれば，べき級数の打ち切り項数S = 15 かつ級数の打ち切り項数J =―33―

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

表-3

第 47 号(平成 22 年 11 月)

外部で液体に接する短い円筒梁 (L / Ro = 2) と長い中実円筒梁 (L / Ro = 40) の振動数パラメータ の収束性に与えるべき級数の打ち切り項数S の影響：H / L = 0.5

L / Ro

S

1st

2nd

3rd

4th

5th

2

5
6
7
8
9
10
11
12
13
14
5
6
7
8
9
10
11
12
13
14

3.50903.50903.50903.50903.50903.50903.50903.50903.50903.50903.48533.48533.48533.48533.48533.48533.48533.48533.48533.485321.55621.55621.55621.55621.55621.55621.55621.55621.55621.55620.40020.39920.39920.39920.39920.39920.39920.39920.39920.39961.04059.65059.65059.63459.63459.63459.63459.63459.63459.63458.00556.69456.68556.67256.67056.67056.67056.67056.67056.670125.36124.99118.20118.20118.01118.01118.01118.01118.01118.01118.92117.79111.59111.46111.30111.27111.27111.27111.27111.27557.98219.02
218.46197.20197.18196.04196.04196.02196.02196.02544.21207.71205.41
186.14185.81184.85184.79184.78184.77184.7740

201 程度まで有効数字十桁以上の積分精度が確保できるRo = 2) と長い中実円筒梁 (L / Ro = 40) の振動数パラメー

ため，本報告ではこれを用いる．なお，構造解析プログラ
ムは，personal computer (Windows XP) の環境下でFortranが示してある．ここで，液高比H / L は 0.5 に設定し，S はによる倍精度計算で開発した．また，全ての数値計算結果

5 から14 まで変化させている．

タ の収束性に与えるべき級数の打ち切り項数S の影響

これより，L / Ro の値に係わらず，S の増大に伴い，基は，有効数字五桁で整理している．
表-1 と表-2 には，それぞれ，外部で液体接する短い中本振動から五次振動までの振動数パラメータ は一定値

実円筒梁 (L / Ro = 2) と長い中実円筒梁 (L / Ro = 40) の基へ収束しており，数値解の収束状態は良好であると判断で本振動から五次振動までの振動数パラメータ の収束性

きる．ここで，本解析法による数値解の妥当性を直接検証に与えるべき級数の打ち切り項数S と級数の打ち切り項

できる参照解は存在しないが，数J の影響と精度比較が示してある．ここで，液高比H / L

1)

は 1 に設定し，S は 5 から 14 まで，J は 11 から 151

本解析法により求めた液体の高さと梁の長さが等しい場合 (H / L = 1) の外部で液体に接する中実円筒梁まで変化させている．また，本解析法の妥当性を検証するの基本振動から五次振動までの振動数パラメータ
は，厳密解と有効数字五桁でほぼ一致すること．

ための参照解が見当たらないため，液体の高さと梁の長さ2)問題の厳密解を求め，これを参照解とすることとした．

本解析法により求めた液体との非接触領域を有する(H / L = 0.5) 外部で液体に接する中実円筒梁の基本振表-1 と表-2 より，L / Ro の値に係わらず，S と J の動から五次振動までの振動数パラメータ の値は，Hが同じ場合の外部で液体に接する中実円筒梁の自由振動/ L = 1 の結果よりも大きくなっていること．

増大に伴い，基本振動から五次振動までの振動数パラメー
タ は一定値への安定した収束状態が示されている．し以上，二点から判断すれば，本解析法による数値解は妥当

かし，中実円筒梁が長いほど，J の値を大きく取る必要がな結果を与えていると言えよう．従って，本解析法によれば，任意のL / Ro とH / L を有する外部で液体に接する中あるようである．ここで，得られた収束値に着目すると，その値は，厳密解 (Exact) の結果と有効数字五桁でほぼ一実円筒梁の固有円振動数の高い解析精度が得られる．また，
L / Ro とH / L を変化させた種々の詳細な検討を行った結致しており，本解析法による収束値の妥当性とその高い解析精度が確認できよう．さらに，2 ≤ L / Ro ≤ 40 の範囲で果，べき級数の打ち切り項数S = 15 と級数の打ち切り項数

あれば，J = 151 を用いることにより，外部で液体に接す

J = 151 を用いれば，2 ≤ L / Ro ≤ 40 の範囲である任意のH /る中実円筒梁の固有円振動数を統一的に解析できる．

L を有する外部で液体に接する中実円筒梁の固有円振動

続いて，表-3 には，液体との非接触領域を有する外部で液体に接する外部で液体に接する短い中実円筒梁 (L /数を統一的に解析することができる．よって，以後の数値計算では，これらを用いる．

―34―

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

表-4

第 47 号(平成 22 年 11 月)

外部で液体に接する中実円筒梁の振動数パラメータ に与える長さ－外径比L / Ro と液高比H / L の影響

L / Ro

H/L

1st

2nd

3rd

4th

5th

2

0
0.5
0.8
1
0
0.5
0.8
1
0
0.5
0.8
1
0
0.5
0.8
1
0
0.5
0.8
1

3.51603.50903.43693.29473.51603.50103.3684
3.14393.51603.49633.33623.08273.51603.49023.30023.02073.51603.48523.27592.981922.03421.556
20.66920.46122.03421.08419.957
19.63122.03420.85019.73119.29322.03420.581
19.54418.92922.03420.39919.45718.68961.69759.63458.80858.49761.69758.10156.56955.90361.69757.51355.53654.69561.69756.96254.43053.35761.69756.67053.77952.483120.90118.01116.94116.14120.90115.51112.87111.27120.90114.23110.68108.65120.90112.57108.00105.48120.90111.27106.14103.29199.86
196.02194.44193.54199.86191.85188.46186.16199.86189.55184.97
181.72199.86186.72180.41
175.91199.86184.77177.00171.59

6

10

20

40

最後に，表-4 には，外部で液体に接する中実円筒梁のる検討を行い，加えて，本解析法を用いる際に現れる幾つ

基本振動から五次振動までの振動数パラメータ に与えかのパラメータをどのように設定すればよいかについてる長さ－外径比L / Ro と液高比H / L の影響が示してある．

も明らかにした．その結果，数値解の収束状態は良好であ

ここで，L / Ro は 2 から 40 まで，H / L は 0 から 1 まで変化させている．自明であるが，H / L = 0 は，液体に接り，また，高い解析精度が得られることが確認できた．さらに，2 ≤ L / Ro ≤ 40 の範囲であれば，べき級数の打ち切していない場合を意味する．

り項数S = 15 と級数の打ち切り項数J = 151 を用いること

これより，外部で液体に接する中実円筒梁の基本振動かで，液体の高さに係わらず，外部で液体に接する中実円筒ら五次振動までの振動数パラメータ の値は，L / Ro に係

梁の固有円振動数を統一的に解析可能であることを明ら

わらず，H / L の増大に伴い減少する．これは，液体の付かにした．

加質量の影響により，振動し易くなることを意味している．
また，L / Ro の値が大きくなり，中実円筒梁が長くなるほ

種々の数値解析法や今後開発される数値解析法の数値解

ど，振動し易く，液体の付加質量の影響を受ける．

のベンチマークテストに使って頂ければ幸いである．

５．あとがき

参考文献

本報告で示した表の値は，既存の有限要素法を含めた

1)

桜井彰雄：水中に立てられた柱状構造物の振動，土木

2)

技術，第16巻，第6号，pp.11-17，昭和36年．
Wang, C.M., Reddy, J.N. and Lee, K.H.: Shear deformable

本報告では，べき級数を試行関数に採用したRitz法(p-Ritz法) を用いて，任意の液体の高さを有する外部で液体に接する中実円筒梁の自由振動問題を定式化し，実用的

beams and plates – Relationship with classical solutions,

な構造解析プログラムを開発した．また，本解析法によりElsevier, 2000.

求められる固有円振動数の収束状態やその妥当性に関す―35―(2010.9.30 受付)