九州大学学術情報リポジトリ

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

一自由度Duffing形自由振動方程式の厳密解の計算法 に関する研究

岡部, 匡

https://doi.org/10.11501/3111012

出版情報：Kyushu University, 1995, 博士（工学）, 論文博士 バージョン：

権利関係：

ウUAUT

第3章 対称Du伍ng系の自由振動の厳密解とその高精度数値化

3.1 対称Du白ng系の自由振動の厳密解 本節においては， 対称Duffing系，

今

_at� +αx+ y X3 = 0，

1

_?

ただし，y �O，一∞<α，y <∞

j

^(3.1.1)

の自由振動の厳密解について包括的に論じることとし，対称Duffing族に含まれる 個々の振動系(後述〉に対する厳密解の高精度計算法については， 次節以降に個 別に詳論することとする.

まず， 式(3.1.1)を無次元化すると次式を得る.

d2X -7 +LXX +CX3=O，周期(T)= 2πIQ

dT" ^{A A} (3.1.2)

ただし，

Lx三sgn(α)， Cx三sgn(y)， X三 x

þ疋1

_，

^ì

日

記

_。三ω1

M ^j

あるいは，

(3.1.3)

Lxニ0， Cx = 1 ， X =

刈

_{y ，}

ì

T= t， Q=ω，1:こだし，α= 0， y ^> 0

j

^(3.1.4)

ここで，ωは(固有〉円振動数を表し，2π/ωはtに関する周期，ω12πは(1司 有〉振動数を表し， いずれも実次元量である. また， 無次量Tに関する周期jは 2πIQである. 実次元量 x， t，ωと無次元量: X，T，Qの対応は式(3.1.3 )， (3.1.4 )

による . 容易に知れ る ように X=玄(T) が式 (3.1.2 ) の解で あ る とき，

X = :1:X(T+ �)， (ζは任意の定数〉もまた解である.

式(3.1.1) あるいは式(3.1.2) において， ばね関数中に現れるパラメータ仏Y (Lx'Cx)の組み合わせにより， し1くつかの代表的な非線形ばね特性が表現され

る. この振動系において周期解が存在するためには，α，y(Lx'Cx)の符号の組 み合わせには制限が生じ， それらを表3.1 に示した. 本論文においては，

1J A『

α> 0^{， r} > 0 ^(Lx = ^{1， C x} = 1 )の場合は対称ハードスプリング系(Symmetrical hard spring system)， α> 0^{， r} < 0 ^{(L x} = 1， ^{C x} = -1 )の場合は対称ソフトスプリング系 CSymmetrical so食spring system)， α< 0^{， r} > ^{0 (Lx二-1， Cx} = 1 )の場合は対称スナ ップスルースプリング系または対称飛び移りばね系(Symmetrical sn ap-through

spring system)，α= 0^{， r} > 0 C ^{L x} = 0^{， C x} = 1 )の場合は単項3次曲線ばねCSymmctrical single term cubic spring system)と呼び分けることとする. 図3.1に はこれらの振動 系のばね関数を示した . 特に対称スナップスルースプリング系については， 初期 条件の違いに より二様の周期解が存在し， それぞれを両振りモードCPu ll swing mode)および片振りモード CHalf swing mode)の解と称することとする. 両振り モードの解は， ばね関数の変位の原点CX = 0)に関し対称となるように振動を 行うものであ り， 一方片振りモードの解とは， 2個存在する安定平衡点のうちの ひとつ(X=1または-1)を変位の中心として 振動 する非零定符号の解である.

図3.1のばね関数の簡単な視察から ， 各振動系の解の存在範囲を次のように特定 することができる.

1)対称ハードスプリング系:O<mplX(T)|<∞

2)対称ソフトスプリング系:O<mfxIX(T)I<l 3)対称スナップスルースプリング系

・対称両振りモード;d<mpx|X(T)|<∞

・対称片振りモード;(〕<mpx|X(T)|<ゾ2

これら 対称Duffing族に含まれる振動系の解は2.1節で論じた3種のヤコビの楕 円関数を用いて記述することができる. まず， 3種のヤコビの精円関数を一指し

て 次のように表示する.

e(u)三snu^U cnu ^U dnu， こ こ で， max(u) e(u) _{= 1 (3.1.5)} 式( 2.2.7) に て示した ヤコビの楕円関数の微分公式を利用すると ， e(u)についての 2階微分は次のように表される.

竺

竺些f11

₌ 一刊(ω仏砂仇卯いl eρe肘+2Ø:>.e3砂

dT耐 1 ^{' .l} (3.1 .6)

表3.1

対称Du伍ng系のばねの分類と名称、

Name of sprmg ^α y

S vmmetrical hard sorin 十 十

Svmmetrical soft sorin +

S

.+

Svmmetrical sin2.1e term cubic sorin2. 。 +

1.5

t-'

\ 1

lγ

Symm巴trical soft spring

0.5 1

1 匂白a

にル-J_-a

4 ⁴ 1

\/

/ーー、‘

" " ，..-'"・・

v 、、

、

-E

^{。 \}

I

/f、\hJ/ グ

^{F 、、\ ノ}

-0.5 _，

ア

-11

Displacement X

，

/

，

、/

， 、1

、1

1

図3.1

対称Duffing系のばね関数

Lx

1 1

-1

。

，

1.5

-44-

Cx

1 -1

1

1

ペJA『

表3.2 ø1， ø

3

_の

定義

e(u)

11

_{ø1 (μ。=0} 山1)

1

ω。=0 ... 1) snu 1+μ。 (= 1... -・.2) 一μ。(= 0 ・・ー1)

cnu 1-2μ。(二 μ。(二 0 ・・・ 1)

dnu μ。-2 (= -2… 1

表3.3 解の適合関数

Spring type

11

^e(u)

1

^μ。

S ymmetrical hard spring cnu o ，_ 1/2 Symmetrical so丘spring snu o

，_

₁

Symmetrical snap-through spring cnu 1/2 ，_ 1

トーーーーーーーーートーーーーーーーーーーーーーーーーー，

dnu ^o ，_ 1 Symmetrical single te口n cubic spring cnu 1/2

上式中に現れるゆl' ø3は， snu， cnu， dnu のそれぞれに対応して定義される楕円関 数の母数 μ。の関数であり， その関数形を表3.2 にまとめた.

式(3.1.2)の解を， 式(3.1.5)の表示式を用いて，

X(T)二Axe[u(T)]， AX(:;t: 0):定数 (3.1.7)

と仮定し， 式(3.1.7)を式(3.1.2)に代入し， 式(3.1.6)の関係式を用いて整却すると 次の関係式を得る.

(Ø1 +2Ø〆).(22二Lx + CXAXγ，

ì

ただし，Axe[u(T)]:;t:O

J

^(3.1.8)

両辺のe に関する 等べき項の係数の等値により， 次の関係式が得られる.

ゆ1 .(22= LX' 2Ø3.(22 = CXAX 2 (3.1.9)

すなわち，

Axz二(2Lx / Cx

)(Ø3

_{/ øJ}

d二Lx /や1

(3.1.10) (3.1.11)

-46-

上式を参照すると 両辺の符号の関係 からLXØ1 ^> 0， CXø3 ^> 0であることが必要と なることがわかる.こ の関係と 表3.1 ， 3.2に示した関係式を参照すると，各スプ

リング系に対応する e(u)の関数形は表3.3に示すように一意に確定することとな り， また母数 μ。 の適合範囲も同時に決定される.すなわち，対称ハードスプリ ング系の解はcn関数，対称ソフトスプリング系はsn関数 ， また対称スナップス ノレースプリング系の解は， cn関数また はdn関数lこより記述されることが判明す

る• e[ u(T)]のTに関する円振動数をρとすれば ，式(2.2.8)"'"'(2.2.10)で示したヤ コビの楕円関数のフーリエ展開公式を参照すると ，

一三一=m: (snu ， cnu)，

KoCμ。) ^{" /'}

I

JL=Z:(dnu)

|

Ko(μ。) ^{， / )}

であればよく，式(3.1.6)のQは次のように定 まる.

(3.1.12)

Q=Qkh(μ。) : (snu ， cnu)

ì

ト( 3.1.13) .a =.QK(μ。)/2 : (dnu)

J

式(3.1.10)， (3.1.11)， (3.1.13)を組み合わせると式(2.2.8)"'"'(2.2.10)からAxe[ u(T)]は 次のように整理される.

Axsnu山品2

24

s^山=山芝�sinnm (3.1.14)

川(白州1^{- q}。 ^凶�d)sinhnQ

Axcnu(T) = 2-fi.a芝

型22L

cosnZ=2JZQZ^{国 l} mnZ(3 1 15)

nE1(odd)1+qon n2ゆdd)coshnQ

Axdnu=J20

1

L Z立

ら

c伽ぽ

1

=-fi.a

l _�十三 l n C刷出

1 2 ^位刷)1 ⁺ qo"'''

I L

2 ^時制)cosh2nQ

I

(3.1.16) さらに ， 式(3.1.13)を， 式(3.1.11)に代入すれば 次のような関係式が得られる.

。2= Lx /ゆl[Ko(μ。)]2 : (snu ， cnu)，

ì

ト (3.1.17)

β2 = Lx / ØJK(μ。)/2 ]2 : (dnu)

J

以上より，式(3.1.10)，(3.1.17)をまとめて厳密解を特定する際の未定定数の決定条 件式を再記すると 次のようになる.

Ax2 ⁼ (2Lx /CX)(Ø3 /ゆJ， (Ølヲム:表3.2参照) e(u ) ⁼ snu ^U cnuの場合; [22 ⁼ Lx /øJ�(μ。)f，

e(u ) ⁼ dnuの場合 ; [22二Lx/ øJKo (μ。)/ 2y

47-

(3.1.18)

上式において， ø1， ø3はともに母数 μ。の関数， Lx' Cxは与えられる定数である.

初期条件として振幅Axが与えられる場合には， 決定すべき未知変数は振動数。

と母数 μ。となり， 振動数。が与えられる場合には， 未知変数は振|幅Axとi手数 μ。となる.

ここで， 式(3.1.15)�(3.1.17)におけるフーリエ係数の計算について若干の考察 を加える. 式(3.1.15)�(3.1.17)には q。およびQの二様の表示を示したが， まず Qからそのフーリエ係数を計算する場合 ， 1/ sinhnQ， 1 / coshnQなる関数計算 の 誤差増幅率を計算すると，

snu : (((1/ sinhnQ) / /Q))二nQ/ tanhnQ →O(nQ)→∞ (Q→∞)

→1ヲ (Q→0)

cnu : (((1/ coshnQ) / /Q)) ⁼ nQtanhnQ →O(nQ)→∞， (Q→∞)

→O(n2Q2)→0， (Q→0) dnu : (((1/ cosh2nQ) / /Q)) ⁼ 2nQtanh2nQ →O(2nQ)→∞， (Q→∞)

→O( 4n2Q 2)→O，(Q→0)

(3.1.19) であり， sn， cn， dn関数のすべてに対して Q→∞ (μ。→0)においては計算精度l二の 問題がおこる. 一方q。を用いてそのフーリエ係数を計算する場合における

ぷ-;;

/ (1- qon) ，

ぷ工

/ (l+qルqon/ (1 + qつなる関数計算に際しての誤差嶋中砕 を計算すると，

一n(1+qJ

S釦nu: (何([旬qon/2/ ( ο1ト一qρ]//勾qo))一 →一コ， (ωqo → ^O_町 ) 2(1一qJ

→∞凡ヲ(ωqo→ 1り)

C∞nuι: (伺([匂qon/2/ (ο1一qon)]/ /勾qo))二n川(1ト一qJ →竺， (ωq2(1 + qon o→O町)

→0， (ω肌q仇O→1り) 一 n(1一q

f

₎

dnu : (([qon / (1- qo2n)] / /qo)) = →n， (q。→0) (1+qo2n)

→0， (q。→1 )

(3.1.20)

48

であり，sn関数に対するフーリエ係数の計算では q。→1において危険となる. ゆ えにこれら楕円関数のフーリエ係数の計算に際して， 計算精度上q。または Q に よる表示式の使い分けを行う必要が生じて くる.その切り替え点は，式(2.1.16) と 同様に μ。=1/2 (qo =ピヘQ =π/2)が最適であり，μ。< 1/2 であ ればqo' ^110> 1/2 の場合にはQによる表示を利用する.

以上より対称 Du伍ng 系の自由振動の厳密解を計算する大筋の流れを述べる と ，まず振幅Ax が与え られる 場 合 に は ， 式 (3.1.18)の第1式か ら陽に

μ。=μ。(Ax)を解く.そして式(3.1.18) の第2式により D=D[μ。(Ax)]を得， さら に式(2.1.9)からQ ， qo = fun.(μ。)を得ると， 式(3.1.14)�(3.1.16)から解X(T)およ

びそのフーリエ係数が決定される.一方ρが与えられる場合には，まず式(3.1.18) の第 2 式から μ。=μ。(D)を試行解法により得 ， 式(3.1.18)の第 1 式からAx' 式 (2.1.9)からQ ， qo = fun.(μ。)を得 ， さらに式(3.1.14)�(3.1.16)から解X(T)を得 る. 図 3.2にはこれらの計算のフローを示 した . 無論これは計算の大筋の流れ の 概略図であって， 厳密解の高精度計算のためには， 前章の議論も含め た詳細!なる 計算式の吟味が必要 となる.

Frequency

Q Q

州ムAAXM

Eqs.(3.1.18)， (2.1.9)， (2.1.13)

μi' qi' Q

Eqs. (3 .1.14) � (3.1.1の Fourier coefficients

Solution ， X (T)

図3.2 対称Duffing系の自由振動の厳密解の計算の流れ

-49-

3.2 対称ハードスプリング系

本節では， 前節までの議論をもとに， 対称ハードスプリング系の自由振動厳統 解の高精度計算法について考察する. この振動系の運動方程式は式(3.1.2)および 表3.1より次式で表され，

-d2X dTτ" +X +X3 = 0，周期(T)= 2π/ο (3.2.1)

その解 X(T)は， cn関数を用いて次のように表記される.

X(T) = Axcnu(T)， Ax = max(T)IX(T)1 ;定数 (3.2.2)

母数 μ。は， 表3.3 より んご0"""1/2の範囲に制限されるため， 2.1 節で論じたよ う なμ0'μ1の使い分けの必要はない . 対称ハードスフリング系の解は， 対称 Duffing族の自由振動解のうちで最も性質が素直であり， 振動は唯一存在する安定 な平衡点(X二0)を中心として実現し， その振幅は0"""∞までのすべての範囲で 存在する. 未定定数の決定条件式は， 式(3.1.8)および表3.2を参照して， 次のよ うに表される.

Axz=2μ。1(1- 2μ。)，

1

ここで， μ。=0"""1/2付Ax = 0"""∞J

。2= 11 Ka\1-2μ。)，

l

ここで， μ。=0"""1/2，件。 =1"""∞J

これらから陽に現 れるんを消去すれば次式を得る.

1+Axz=ζ2�2，

1

1+AX2~02，ただし， (Ax →0， �→ 1) j

3.2.1準線形微小振動(Ax →O， �→0，μ。→0) 1) 振幅Ax(→0)が与えられる場合

(3.2.3) (3.2.4)

(3.2.5)

振幅Ax が初期条件として与えられる場合には、式(3.2.3) を変形した次式

μ。=Ax2/2(1+Ax2) (3.2.6)

から、μ。(Ax )を得る. これを2.1節で示した式(2.1.11) または式(2.4.6)に適用し て，

ハU CJ

qo = qo[μ。(Ax)] (3.2.7)

を 得， その途中で定まる K。と入力値Axを， 式(3.2.4) に代入して0を求める.

さらこれらq。と ρを式(3.1.15) に適用すると， フーリエ係数が計算され最終的 に解X(T)= Axcnu(T)が決定される. ここで， ο→1における.0= .o(qo)の高精 度表示として (.0-1)を得 るには， 式(3.2.4)を変形した次式を用いる.

(.0 -1) 1 [1 + (.0 -1)]] = 1-11β

=2 μ。-4Pa(1 + 2Pa)(1-2μ。) 1+κ

�

^1-2μ。

= 0(12qo)' (q。→0) (3.2.8)

すなわち，

β-1 = (1-11ρ) 1 [1-(1-11ρ)] (3.2.9)

ここで， 最終的に厳密解が特定されるまでに用いられる一連の計算 式に対して，

Ax →0の場合についての誤差の増幅率(2.5節〉を計算しておくと，

式(3.2.6) ; ((μ。11Ax>> = 2/(1+Ax2)→2， (Ax→0) (3.2.10) 式(2.1.1 7) ;

((q 01 1μ。>>→O(D，

1

・L.>.

ト

^(3.2.11)

((Ko Ilqo>>→O( 4qo)→0，[式(2.5.4)，(2.5.5)， 図2.7参照]

J

バ(3.2.4) ;

((.01 1Ax>> = Ax2 1 (1+Ax2)-((民11μ。>>/((Axllμ。>>

→0(1)， (Ax→0)

、...

^{，，_・、戸}

1.- 1.- Y'-， (3.2.12)

((Ko 1 1μ。>>→O(4qo)→0， [式(25.4)，(255)，図(2.7)参照]，

((Ax 1 1μ。>>→0(1) 式(3.2.8)，(3.2.9) ;

(((.0-1) 1 Iqo>> = (((.0-1) 11(1-1/.0)>>(((1-11β)llqo>>>

→1， (Ax→0)

、Elll、pti--1ノ 司コ〆'E、、 っ“ 114 、、，/内4〕

ここで，

-圃a・

ウ山戸、J

であるので， θ→1に対しては計算精度上の問題はない.

3.2.2 局周波振動(Ax→∞，f2→∞， μ。→1/2) 1)振幅Ax(→∞)が与えられる場合

この場合低周波振動の場合と同様に， 誤差の増幅に関しては安全であるのでιに (3.2.3)， (2.1.17) [式(2.4.6) ]を用いてμo，qo'�。を計算すればよい. ただし，式(2.1.17)

より求まったq。を用いて， 式(3.2.8)からβ=f2( qo)を計算する場合 ，

<<(f2-1) / /qo>> = <<(f2-1) / /(1-1/ f2)>> <<(1-1/ β) / /qo>>

→∞， (Ax→∞)

"... "... 1、...

」ー'--V'-， (3.2.17)

< < (1-1/ .0) / / q 0> > →∞

<<f2-1/ /(1-1/ .0)>→∞

となり， Ax →∞では危険となる . ゆえにAx→∞の場合に0を求めるには， 式 (3.2.8) によらず式(3.2.5) によって， f2=θ( qo)を計算する . 式(3.2.5) の誤差の増 幅率を計算すると， <<β/ /Ax>>→0(1)， ( Ax→∞)であり， Ax →∞では安全である ことが確かめられる .

2)振動数(.0→∞)が与えられる場合

。→∞の場合には， ρ→1の場合と同様の手順にて計算を実行すればよいが，

求められたんから， 最終的に振幅Ax二Ax[μ。(.0)]を式(3.2.3) を用いて計算する 場合には

<<Ax / /μ。>>=1/2(1-2μ。)→∞， ( μ。→1/2， f2→∞) (3.2.18) となり， 。→∞に対しては危険である. その ため， .0→∞に対しては， 式(3.2.5) を用いてAx= Ax(f2)を計算する. そのときの誤差の増幅に関しては

<<ρ/ / Ax)) = AX ² / (1 + AX ²) - < <KO / /μ。))/ <<AX / /μ。))

→1， (.0→∞) (3.2.19)

となり安全となる.

以上より， AX，f2 いずれが与えられる問題に対しても ， ともかくAx，f2→∞ な

ーーー圃・圃圃・

-53一

る場合に限り， 式(3.2.5)にKo(qo)，Ax'nを用い， Ax，nを得るのが一つの留志点 である.

3.2.3 対称単項3次曲線ばね系との対応

式(3.2.1)において， 線形項を欠く場合〈α=0 ，対称単項3 次曲線ばね系〉に 対しても， その解は 式(3.2.2)にて記述することができる. ただし， その場合， -J手 数等のパラメ?夕は次の値に固定される.

μ。= 1/2， qo = eπ，Q=π/2 (3.2.20)

そのとき， 振幅Axと振動数。の関係は次式にて表される.

Ax = Ko(1/2)Q = 1.180 355 759 04…xQ (3.2.21)

ここで， 対称ハードスプリング系のスケルトンカーブは0→∞において， 対称 単項3次曲線ばね系のそれに漸近するが， その対応を考察しておく. ここで， 適 当な誤差評価基準による有限のAxあるいは0より大なる場合をα=0に対す る特性と見なす. まず， μ。三1/2 の場合の相対誤差をCm(> 0) とおけば， 次式カ、

成立する.

μ。=

j

^(1ーら)

上式と式(3.2.3)， (3.2.4)から ， Axz二 (1-Cm) /εm信1/cm'

.02 =1/凡[(1/2)(1一九)]εm ただしん→0 ゆえに，

Ac三1/

fi:

， .0二Qc三1/�(1/吟

Fm

(3.2.22)

を得る Ax>>Ac， Q>>Qc なる直線 部分は ， μ。の相対誤差がCm であるときの 対称単項3次曲線ばね系のスケルトンカーブと考えることができる. すなわち，

μ。の相対誤差がιならば， この部分を対称単項3次曲線ばね系と考えてよい.

さらに μ。の相対誤差がら であるときの qo(μ。三1/2) の相対誤差を考えると， 式 (2.1.10)， (2.1.11)より，

μo�l6qo/�\ Ko信(1⁺ 2qO)2， ただし， μ。→1/2

-・圃・

AUT CJ

の近似が成立し，

<<μ。//qo))店1-^&ZO / .J Ko， ただし，μ。→1/2 であることより， q。 の相対誤差は，

ん/[1-8e-tr /

，;

K(l/ 2)] � 1.46677 2…Emヲ ただし，μ。→1/2 (3.2.23)

となり， この際のフーリエ係数[式(3.1.15日の相対誤差はιにより評価される.

以上のように， Ax >>Ac ， .(2)>.(2cにおいては， 基礎式(3.2.1)のパタメータαの 変域は事実上零まで含めることができ，α=0とα:;t:0の場合を統一することがで きる.

3.2.4 フーリエ級数の打ち切り誤差

式(3.1.14)�式(3.1.16)にて厳密解を計算する場合， このフーリエ展開の無限級 数和は実際の計算においては有限項で打ち切らざるを得なし1から， それに伴う誤 差を検討しておくことが必要となる. とくに後述するようにソフト， スナップス ルースプリングスプリング系の場合の超低周波領域において高精度な近似解を作 成するためには， 必要となるフーリエ級数の項数も膨大な数に達するため， その 項数をあらかじめ把握できれば， 実際の計算において好都合である. 本的jでは，

式(3.1.15)のcn関数に関するフーリエ展開の有限項打ち切りによる近似誤廷の評 価法について検討を試みる.

まず， 議論の出発点として， 次なる関数を定義する.

Pc

三三n孟5Ed的ω)

δ D ι n 三 すγ 回

₁

一 n?品芯:石dod ) C∞05油hnQ )

(3.2.24)

ここで， 無限級数で定義される厳密解[式(3.2.2)， (3.1.15)Jを次式で示すように有 限級数にて近似した場合を考える.

川 一 叩

d 。 けでム問

g d一川

d 。 JZ附

T X Tt x

_(3.2.25)

有限項打ち切りによる近似解の絶対誤差は，

-

戸、J戸、J

|防δ汀Xん(仰川N

相対誤差εは，

(3.2.26)

E二|ぽ(N)(T_)

I = IX(N/T) - X(T 2 1

^<

1 δPc

max(T) X(T) I ^I max(T) X(T) ^{I I} Pc

^(3.2.27)

で表される. ここで， 次の近似が成立する.

|δX(N) (T)I = 1

^2J2Q

� _{1 ど cosng - Y cosnj 引}

1

^n(ぱ;1)

co

^sh

nQ

^n�)

cosh nQ ^J

= ¹

²

^J2.ρ〆N

_{n川(ぱ叫d)}

cosnn

^{ια r}

p c

r AU σ Q っ“ 一 Q

一n1i一hu一nu一戸lv

O L寸ムM

ウIH /削V ウム <一

_(3.2.28)

ImaX(T)X州三ldo Z) _COS L _Q|=ld今cl

^(3.2.29)

式(3.2.24)のPC' oPc の近似値を得るため，

Euler-Maclaurinの公式(138)を用いる Euler -Maclaurinの公式とは，

(無限〉和を積分近似により簡単な数式で表現しよう とするものであり， 次式にて表される.

土f(α+幼)三�r

_M付^+nh

f (x)ゐ+11f(α)+ f(α+nh)l

2L- ^{， / - ， /J} 式(3.2.30)を用いると， 今c の近似値は，

(3.2.30)

δ Pc こ

_n

ど

_�)

coshnQ

¹

^{= Z}

^�

cos

^h

^{[ (2k + N +} 2)Q]

1 r∞ dχ 11 1

三 一一一 I 一一一一+一1 ^----

2Q 削+2) Qco

^sh

^x

²¹

cos

^h

^[

⁽

^{N +}

^2)Q]

1

1

+ 1 t釦叱e 一( 川N仙川 + 2叩2

2c∞osh[(N

+ 2勾)Q] Q

^{， / 1}

+eピ一4引2(げN+リ2)Q

必; cxp[(N+

_{1 +}

_{exp[2(N +}

2) ₂₎

ln ln ^-Jqo

_V

_{- -L} .Jqo] _u/-， ^] ₊

^.

_t

^.

an

^--}

^{1 {exp[(N +}

²⁾

ln ^..Jqo]}

一

_ (N+2)/2

= 一 l

^In1

^J

^{v q} � ⁰ v 1 ^{q'1.0 +}

q

⁰

川伽

�" ."

+ 斗2 ")+

，

_t an

^正

^ぺ

^-

' ^{q 仇} �

^v

^{0f 2 ) 川 勺 / 2 2 )}

^(3.2.31)

同様に，

一-

〆り宍J

Pc=1-i-+1tan-l(e-Q)二f-QZ+1t加\e Q) 2 coshQ ^Q

r exp[ln刊。]

ゾqO V 1U +tan正-l[ exp(lnゾqo)月] ゾ

刊仙q仇o

l+e叫X勾p[仰2引ln

必 ; h ^]

二 一叫

市;

^{+tan-1(Jqo) (3.2.32)}

式(3.2.31)，(3.2.32)より， 式(3.2.27)のεは次式にて表される.

� I� …η + lぱ (ピ(N+2)Q) 初，. 1 11+e-��/VH)k! ⁰

ε=I� I=I � ~

Pc

I I

一三二一+1tan-1(e-Q)

1 +e-2Q ^Q 一(N+2)!2

-ln

J

qo qo _仙2)+t^an

吋

q0 (N+2)β)

v au1+qJ J u

(3.2.33) ー叫

市;

^+tan-ld)

准線形振動域(q。→O，Q→∞，.(2→1，Ax →1)では，

Pc三1/2coshQ→〆(=^�qo)'

δPc ⁼⁼ 1/2c

い 恥

e-^(N+2)

刊

lⁿ)^，

}

^(3.2.34)

E三|δPc/ Pc I→e-(仙刈=q。(N+1)/Z

であり， 要求される相対誤差限界εを満足するために必要となるフーリエ級数の 項数Npは，

N ^� -ln c / ^{Q - 1} = ln c / ln長-1 →Np三lnc/ ln長+1 _ì

Nミlnc/ In(A / 412)-1 →Np三lnc/ In(A / 4-J2) + ¹

�

^(3.2.35)

N三lnc/ ln(必弓/213)-1 →Np三lnc/ In(.J万�

1:

_{/213) + 1} ^j となる.

また高周波振動域(Q→π/2，ο→∞，Ax →∞)では，

Pc三eπ匂1+2/Jr)， _ì δPc三e刊均吋1+2/π)，ト

ε=1δPC / PC I→e -(N+l)rr/2 j (3.2.36) であり， 相対誤差限界εを満足するために，

一-

勺J戸、J

N三-2lnE-1 →NF三-21nE+1

π π (3.2.37)

なる項数を必要とする.

図3.2には， E = 10-5， 10-14， 10引を満足するに必要なフーリエ級数の引数NF の計算結果を示した 高周波振動域では， 各εに対し，

NF三9 (ε= 10-5 '"'-'単精度丸め誤差) ì

NF三22 (E = 10-14 '"'-'倍精度丸め誤差〉

ト

Npさ48 (ε= 10-32 '"'-' 4倍精度丸め誤差) j

なる項数が必要となる. 一方準線形振動域では， NFは1に漸近していく.

(3.2.38)

100

ε=10-32

�

「

8 g ⁶⁴

I I

^10-14

・比口らaロロ4 ^B 4

10

l=-I

^10-5

匂0⁴

てδ

^コ

1

1 10

Frequency Q

100

図3.3 フーリエ級数の打ち切り誤差εを満足するに必要な項数NF (対称ハードスプリング系〉

一一-

。。戸、J

3.2.5 数値計算結果

図3.4 には， Q ⁼ 1.01， 1.1， 15， 2.0， 25 の各振動数における振動波形をぶした これら の振動波形は0→1では単弦波に 近づき ， 振動数 の増加につれて高調波合 有率一定の歪んだ 波形 となるのが特徴である.

図3.5に 振動数/振幅 の周波 数特性(し1わゆるスケルトンカーブ〉および調波 振幅の周波数特性を示した . 調波振幅 Aの定義は次式による.

， 。。

X(T)二日n^COSF1Z (3.2.39)

n:?:l(odd)

対称ハードスプリング系のスケルトンカーブは， 古典的非線 形振動論による近似 解法によれば，

1 + (3/ 4)Ax2 二Q2，(ただし， 0 < Ax <∞，1<β<∞) (3.2.40) で表され， 周知の簡潔な双曲線となる. 図3.5中の細破線は式(3.2.1)において α=0である対称単項3 次曲線ばねのスケルトンカーブを示し，式(3.2.21)より科 易に 知 れるようにこの系のスケルトンカーブは原点を通る直線で表され， 対杉J\ノ\

ードスプリング系のスケルトンカーブは ο→∞に てこの直線に漸近していく.

各次数の調波振幅は， 。→0では，

Af

dO

_→4d

l tr

_，(ρ→0)

cosnnu ^{\. 11， ノ} (3.2.41)

でなる漸近特性を有する. 一方ο→∞においては，

2-J2Q → 2-J2Q ，(ο→∞)

coshnQ coshnJr /2 (3.2.42)

であり， さらにフーリエ級数の次数n が大なるときには， 次なる近似式が成立す る.

An→4.fiQ(e叩)"，(Q→∞， n，十分に大きな整数) (3.2.43) このため， ο→∞における図示の対数表示では， すべて の次数の振幅特性が平行 な勾配1の直線群となるのが特徴的である. 前述したように 対称ハードスプリン グ系の振動波形は， 。→1では単弦波に漸近し， また 0→∞では歪率一定の波 形に漸近し振幅は振動数に比例して増大すること となる.

-59-

図3.6には母数ノミラメータμ。の関数として， 各種の特性値をぶした. 振r�Mおよ び振動数の漸近特性を母数μ。の関数として示せば， 式(3.2.3)， (3.2.4)より，

Ax →∞，.(2→∞， μ。→1/2では，

4→

作(

^1/

;

_-μ。

₎ I

九川J1

^/2一μ。

^{) J}

^(3.2.44)

であり， 一方Ax→0，.(2→1，μ。→0では，

Ax →

兵五0

^'

¹

.(2-1→3μ。/4 J なる漸近特性を有する.

(3.2.45)

3

� ¹

肖

� ⁰

0

ち -1コ c/)

-2

司、〕

0 π 2π

Independent variable m

図3.4 対称ハードスプリング系の解の波形

//

10-66- / 10-76-' 10-8長一 /

-60-

。。

�/

xσ)=ヱAncosnT

n�l

，r曲q^、、

ノ\

⁾

d

ω 民事弓ご8

、\

ノミ

t r・画、

/

ノ U 与 ³

/ / �/

/ / γ

I / .rI1y 寸/ ぷ// f

I I / .<.J.1シr 寸 〈ザ'/ /

1 100

_。

1

₂ 3

Frequency n ^-

1

Frequency n

図3.5

対称ハードスプリング系の振幅の周波数応答

-61-

0.001

1/2-μ。

0.01 0.5 0.1

50 Eーナナ

10

1

0.1

0.01

0.001

sdJlqhq

hH布、

10

0.5

。 ^ープ/

Ax "〆〆〆 /

/ 記;イふ /

〆イム 〆 _/

。ーシ'f /�.

^d片、

、dG/

』ヘデ

句人/ V�:

イ: > #

1

0.1

0.01

0.001

0.001 0.01 0.1

Modulus parameter μ。

。貯ぷuJlqhqhH可

母数μ。と諸特性値の対応

ヌ13.6

-62-

3.3 対称ソフトスプリング系

式(3.1.2)および表3.2より対称ソフトスプリンク守系の運動方耗式は次式で・ぶさ れる.

d 2X

-dT" 7+X- X3=O， 周期 (T)ニ2π/n (3.3.1)

この振動系には1個の安定平衡点( X= 0)と2個の不安定平衡点( X = :1:1)が 存在し，安定平衡点を中心とした振動が実現する. 解は IX(T)I三1なる範囲のみで

存在し，振幅がIAxl>lとなる場合には周期解は存在しない. 表3.3で示したよう に，この振動系の解は次式のように sn関数で記述できる.

X(T) = Axsnu(T)， Ax二max(T)IX(T)I;定数 ，

ì

Ax = 0'"'-'1件。=1'"'-'0 j (3.3.2)

この ときAx=0~1に対応する母数の 範囲は0<μ。<1であり，2.1節で論じたよ うに一連の計算の定式化に際してはμ。 とμ1の使い分けが必要となる.

まず，式(3.1.18 )から次の 関係式を得る.

AX2=2μ。(1+μ。)，

ì

ここで，μ。=0'"'-'1付Ax⁼ 0'"'-'1 J

。2= 1/ (1 +μ。)K02，

ì

ここで，μ。=0'"'-'1件。=1'"'-'0 J

あるいは，式(3.3.3)，(3.3.4) よりんを消去すると次式を得る.

AX2/2+K0202=1

特に μ。→o(準線形微少振動〉では，

Ax 2/2+β2'"'-'1，

ì

ただし，Ax →o，n→1，μ。→o

j

なる関係式が成立する.

3.3.1 準線形微少振動(Ax→0， n→1， μ。→0) 1)振幅Ax(→0) が与えられる場合

(3.3.3) (3.3.4)

(3.3.5)

(3.3.6 )

この 場合 ，0<μ。<1/2に対して，その振幅は0<Ax <.J2 / 3ニ0.816497・…が対応

「，3fo

表3.5 式(3.3.13) の係数表

λ= 1.416856 _! 12 = 1.728053 _! 13二2.039555

.一一一一一一一一一一一一一一ー一一一一一一一一一一一一一一ー^一 λ= 1.241502 _! 15 = 3.727125 _!

する. 式(3.3.3) より次式を得る.

μ。=Ax2j (2-AX2) (3.3.7)

これを式(2.1.17) に適用し試行解[または式(2.4.6) を利用]により対Jぶする ζ，q。を得 ， 与えられたAxと求まった民，q。により式(3.3.4)， (3.3.5) からθが 計算される. 。 →1 であるから， この場合に高精度の値を得るには次式を川い る.

-[2 (二1-κ

片

⁾

8�(1+2�)+μ。Ko2 8�(1+2�)+Ax2 /2 -

Ko j1弓�[l+Ko j1弓;]-K。(κ+�1-Ax2

/2)

=o

〈

^μ。)= O(向。)， ただし， ([2→1) _(3.3.8) さらに， [2，q。を式(3.1.14) に用いればフーリエ係数が特定される. これら -l!_liの 計算式に対して， Ax →O(μ。→0)における誤差増幅率を計算すると，

式(3.3.7) : ((μ。//Ax>>二4/(2 -AX 2) = 2(1 +μ。)二0(2) 式(2.1.17) : ((KO / /qO)) = 式(2.5.4) →0(4qo)→0

((qo / jμ。)) =式(2.5.4)→0(1)，ただし， (qo'μ。→0) 式(3.3.8) : (((1- [2) / /μ。))=_{V "} � �+(1+μ。)(

α

^。

川

1-

[2 2(1 +μ。):lli.Kol.

(3.3.9) (3.3.10) (3.3.11)

→0(2/3)， (μ。→0) であるからいずれも誤差の増殖に関して は問題ない.

(3.3.12)

2)振動数β(→1)が与えられる場合

0<μ。<1/2に対して ， 0.692…<ρ<1が対応する.振動数。が与えられた場合 には， 直接βを用い るのではなく， 1-[2が与えられる問題と考え， 式(3.3.8) を

-64-

試行的に解き， 対応するKo，q。を得る. この場合の試行解の初期値としては， 式 (3.3.8)の数値的逆関数を最小自乗誤差処理して得た次の多項式を利用すること

ができる.

1 _ r J 5

qo = 二二三[1+:Lfn(l-.oy] _(3.3.13)

12 ^{L ;;;.}

上式の多項式の係_{数えを表 3.5にまとめた}• .o，KO.，q。を式(3.3.3)に適川してAx が求られる . これらのアルゴリズムについても誤差の増殖に関しては式(3.3.9)'"'"'

(3.3.12)により問題はない.

3.3.2 強非線形低周波振動(Ax →1， .0→0， μl→0) 1)振幅(Ax→1)が与えられる場合

Ax →1においては， 対応する母数はμ。→1となるので， 2.1節で論じたように

一連の計算_{に 対しμ。， q。基調によるのではなく， μ}l' ql基調に アルゴリズム を変

更する 必要がある . 0<μ1 <1/2 (・.・1/2<μ。<1)に対しては， -J3/2<Ax<1が刈­

応する . この場合Ax そのもの が与えられるのではなく(1- Ax) が与えられる も のと考えて， 式(3.3.3)から次式を得る.

μ1 = 2(1-A x 2) / [1 + (1-A x 2)] ì _{ト (3.3.14)} ただし， 1-Ax2 = (1-AX) [2-(1-AX)] j

このμ1 を式 (2.1.17)に適用 し試行法[または 式 (2.4.6)を 利用]に よ り，

q1'民= fun.[μ1(1-AX)]を得， 式(2.1.13)よりQを得る . 他方， 式(3.3.4)，(2.1.9)カ ら導かれる 次式，

β= 2Q/πK1

fiヲ;

^{(3.3.15 )}

により0を得る . ある い は式(3.3.5 )か ら導かれる

。= 2Q

�

^{1-Ax2 /2/ 1CK1 (3.3.16)}

を利用しでもよい. 以上より0ニ.0(1-Ax) が特定され， よって(.o，Q)の関数と して， フーリエ係数が特定される. フーリエ係数を計算する場合， q1 では_な く _Q による表示を用いて計算_{することに 注意する .}

cJ fb

以上の 計算においての誤差の増殖に関する危険度を評価すると，

式(3.3.14); (((1-

Ax

2) / (1-

Ax

₎₎₎ = 2/ (1 +

Ax)

= 0(1)

((μ1/ (1-

Ax

)))ニ((μ1/(1-

Ax

2)))(((1_

AX

2) / (1-

Ax

)))

三[1/(2 -

Ax

2)][2/ (1 +

Ax)]

= 0(1) (3.3.18) さらに， 2.5節の式(2.5.4)， (2.5.5)より， q1'

�

= fun.(l-

Ax

)の計算までは安全で^あ

、

(3.3.1 7)

り， また， Q= Q( qJ についても， 式(2.5.7) から， q1→0では誤差の増別に闘し ては問題ない . ここで ，式(3.3.15)を利用して Q=Q( qJを計算するときの誤差増 幅率を計算しておくと，

((β/ /q)) = ((Q / /q)) -

((K，

/ q)) + μ_1 ，((μ1/ /q))

2(2-μ1)

ト

→O(((Q / /q)))→0， ( ql→0) ， [式(25.7)参照]

j

であり安全であるこ とが確かめられる.

(3.3.19)

2)振動数 β(→0)が与えられる場合

0<μ1 < 1/2に対して， 0<β< 0.692…が対応する. 式(3.3.15) から次式をねる.

πQ 1

1 4 →(1)， (μ1→0) 12Q一

五�

)

^{1-μ1/2 (3.3.20)}

Q，μ1 ともに q1の関数であるから， この関係式を， q1 = q1(Q)を試行的に解く問 題とみる. そのために，

L1s三1-JrQ / 12Q，

ì

.・. π2/2Q二(π/12，ρ)(1-L1s)

J

を定義すると，

(3.3.21)

…p( _Jr2 /2Q)二qfls/

い }

ただし， qfls三exp(-π/2Q) (3.3.22)

により， 式(3.3.2 0)は右辺のK"μlとも ql二ql(L1S)の関数， 右辺は1-Asの関数 であるから， 結局次のように変形される.

L1 s (= 1 -1/

K，

�1 - μ1/2)

-66-

)る定 お め 特 ふ 求

も

σ を J

s A

一

と

4

ト

r t ト 一 J

何

時

加 川町 一

軒

芯

1

パDa一

的

淀 レ+一

一 j 匂 用 } 'D 一 /ム ふん 川 リ

μ

一 5 振 一一 一脆引 J 刊 一 ど

、)

山 れ

仲 一 );! 町鵡 A 伎 の ぬ 一

神

清

川 4M 山 一 肌 峨 九町

-

I

1 一2j

川

){一例ワ土勾七 Jdノ

か

用 ら

ゆ ι

一

ゾ リ 2

一

= 位 内 / 適 う + 1 一 寸恥 旬 一 寸 川 J J山 一にM

FF

q 二

+一71 t 2 4:・­づIJ V 一

一円

・

i 2

1 P

一

2

1ノオ)

v h

l M E L t

φ

I H ( A S

L か

τk引fO一2/1hμ

t 一+V九一+H1)+ JLi: 勾 辺 ト 引 符2一1K11てJ一一 引 数 什Iゆ

灯ロ

ベ 式

、モμ一一州関οzc- 一zベ

ぬ

を μ「同L人L一1必 y円 W

ω 内川れι巾+一J

ししm

b吋 一φtし三I

γ… ρ

こ

ρ 土 q一日見一一一。こφ

3 い 'μな0トK 昨 一一一1

f

ベ

これです

た

め

Q

さ 程

で

公

定 過 る パγ

特 中 れ 式

が 途 さ

1-A.. = μ1/2

1-μ1/2十

五ヲ

^1/2)(1-μJ

あるいは，

(3.3.23)

(3.3.24)

(3.3.26)

トAx二_1+、/1-

三

^τ_U ⁼ O(μ1/2)二O(匂.. . ' . -.' 1)

I

� ただし，U三8�(1 + 2�) -L1s (2 -I1s )K¹

2 j

以上の計算においてもμ1，q1→0ならば誤差の増殖については問題ないことが碓 (3.3.27)

かめられる.

さて， 式(3.3.25)より， ο→0(:.Q→O，ql→0)では次の近似式が成立する.

Q→πQ/J2， ・.ql→e-tr/ゐ (3.3.28)

このQ= Q(.o)の近似は， .0< 0.025程度に対して相対誤差 10-75であることが， 式 (3.3.21)， (3.3.22)から確かめられ， また，

Ql""'_'

ムコ

12，...._，μ1/16< 9 x 10てL1s <10-75

-67-

が成立する. このL1sが微小となる領域においては， 式(3.3.23)， (3.3.24)に対し試 行解法を適用したのではtL1s自身の値がアンダーフローをおこすために解を特定 し得ないという困難が生じる. しかし 試行解によらずとも， 式(3.3.26)の近似 式によれば直ちにQ= Q (β)が得られることとなり， フーリエ係数を得るだけな らば問題はない. しかし， この場合は，μl，qI，K1の値を得ること ができないので，

式 (3.3.26)ヲ(3.3.27)の計算に必要なqlを式(2.1.13)を利用して求めることとすれ ば， ο→0 にeおいて危険となる. この場合のql= ql(f2)の高精度計算法を確立す ることは今後の課題となる が， ここでは 式(2.1.13)を利用する際の危険度を把保 しておくと，

((ql//f2))=π/ �2f2 = -lnql

(3.3.29)

式 (3.3.28) が十分成立する f2<O.025 お よびqlのアンダーフロー限界〈ここでは ql ;:::: 10-75と仮定〉に対する危険度[f2 からql(f2)に対する相対誤差増幅率]は，

((ql / /β>) = 88.9""'171 ，

ì

ト (3.3.30)

ただしt ql = 9 ^x 10引""'10てβ=0.025""' 0.013 ， Q = 0.06""' 0.03

J

となり， 単精度計算以外ならば， 実用上かろうじて許容できる程度であることか わかる 上記の β の下限以下(f2< 0.013)の計算は β=0 と同等に扱われてしま うが， 式(3.3.28)なるqlの近似は十分な精度で成立する.

3.3.3 フーリエ級数の打ち切り誤差

対称ハードスプリング系の場合と同様に式(3.1.14)におけるフーリエ級数のイ了 限項での打ち切り誤差を検討しておく. ここで， 仮定として， e-c川)Q< 0.3 [与式 (3.1.4)の フ ー リ エ 級 数 の 打 ち 切 り相対 誤差]と し て お く . この仮定は ， 1/ ^si凶(N+ 2)Q""'2e-CN+2)Qなる近似処理が相対誤差(r._， e-2�N�2)Q )<10%で成虫する

ことを意味する.

式(3.2.24)と同様に次の関数を定義する

N+2 c';�"，rYr N+2 ^c';�...nr

δPS三L �……=2 2;

n

_HG1

L

nCω�) ^sinhnQ nê凶d)e u -e ^u (3.3.31)

-68-

まず， 式(3.3.31)を複素数e'rúffの形で表現し(これを OPsとする)， 式(3.2.30)で 示したEuler-Maclaurinの公式を適用して， 無限和の積分近似によりδPSの近似式 を作成すると，

W.._.，...� _inilI

5九 二2_n(凶d)

乞

_{e -}ⁿ

f

_-e ^{-"η .._}

N+2，…，∞ ，

三2 ) �いν e(j川zω-Q)n + eザeUilIμ仰l河zι山一寸叫3

n 川(0ωdd的)

日(向山)+ん

か

(jilI-Q)x+e(…

n(jilI-QXN+2) n(jilI-3QXN+2) ...，(j町-5QXN+2)

二e(jilI-QXN+2)+し +し +・・・+し +・・・ (3.3.32)

Q- jm 3Q-μ2T 5Q- jm

仮定により， e-(川)Qくく1であるので， δPSは近似的に次のように表される.

δ九 三日ぜeピ汁一

\. Q一〆3万Tノ (3.3.33)

ここで， δPSはδPSの虚部をとることにより， 次のように表される.

δPS三Im(O

Ps)

= e-(N+2)Q

X

{トド討i泊in(頃山(

Q2+べ(m)2川L �----，-

-/--- --- ---，- . • -/--- J

j

= e-(N+2)Q

X

{

^{Sin(N +2ル�Q2+(ぽ)2 1}

^{問 L � ���}

^{l ， �. +， 2-r}

^…

^{-_ ， � � �� n}

^{空 | ト}

^Q

_{J I}

^σ3.34)

式(3.1.4)をn=Nまでの有限項で近似した際の打ち切り誤差δ(N)Xは，式(3.3.34) を用いて以下のように近似計算される.

δ(N)X(T) = X(N)(T) - X(T)

= 2J2.a_I

1

_n(odd) ¹

ど

^Slnn_Slnnn{!

ゲ

^-_n(odd) ^芝

…ゲ 1

_Slnnn{! I ^{= -2ゾ2f2ops}

=

一山 eピ-べ刊小叩(仰仰山叶N仙川+ワ?

| 可ρTII|

+ 1_'> ^{-， __，，>} I sin{(N +2)ぽ+ tan-1

-Q I r

^(3.3.35)

ゆえに，

|δ(N)X I

<212n円切(1+ 1/ Q)

A x= max(T)

X

(T) に対する相対誤差の推定値εは次式により求めら れる.

E二|

δ…X

ームニムー

A X

式(3.3.6)，(3.3.7)よりεを満足するに必要なフーリエ級数 の 項数NFは，

け

^ln

[57 ^(γJ l

^-2

となる. さら に振幅Axは式(3.3.3)から近似的に，

準線形振動; Ax →J32q。→4-/2 e-Q， (Q→∞)

強非線形振動; Ax →1- 4Qo �1- 4exp(ーが/2Q)， (Q→0) であるから， 準線形振動(Q→∞)では，

-ó9-

(3.3.36)

(3.3.37)

(3.3.38)

(3.3.39)

ε →exp[ -(N ⁺ 1)Q] / 2， (n→1， Ax →0) (3.3.40) となる. ここで， Nニ1 とし た場合は， ソフトスプリング系向山振動の基本波の みを用いた単項近似解に相当するが， その場合 の準線形振動域の相対誤差は， 式

(3.3.40)， (2.1.19)を用いることにより，

ε �AX2 / 64� (1- n) /24ヲただし， n→1，Ax →o (3.3.41) であり， 基本波のみを用いた調波ノ〈ランス法， 平均法等の誤差 の推定が容易にnI 能となる.指定された相対誤差εを満足するために必要な項数Npは， 式(3.3.39)，

(3.3.40)， (2.1.19)を用いることにより，

N ^F > ln(2s) /ln(Ax /412) + 1， (Ax→0) ì

(ο→1) J

により求めることが できる.

(3.3.42)

強非線形振動(n→O，Ax→1)では， 打ち切り誤差εは， 式(3.3.37)，(3.3.39)よ り，

ε→(4/π) ^exp[ -(N ⁺ 2)Q] (3.3.43)

であり， このεを満足するに必要な項数は式(3.3.19)，(3.3.28)， (3.3.40)から，

-70-

2

. 8 . 4 .. ")

N F > - 7 in -- ln - - 2 3 (A x →1) I

7rL.

1-Ax

πε

|

_(3.3.44)

nu → ο

ウU4-m

n

ユ一ρ

>

N

となる. 上式は，

振動数×高調波の次数=

O( l ^lnε 1 ⁾

なる関係を示し， 超低周波振動(.(2→0) の振動解の表示には多 くの次数の高調波 を考慮することが必要となることがわかる. 図3.7には， 対称ハードスプリング

系の場合と同様に， ε=10-5jo-14jo-32に対するNFの計算結果を示した

�

.U8 ^ロo ヨ

-ロH4A 102

位Cロ⁴ 匂04

gaB

^ω

100

ο

vJ PしV I n

o∞ GA e T且EA ¹

図3.7 フーリエ級数の有限項打ち切り誤差εを満足するに必要な項数NF

〈対称ソフトスプリング系〉

- ⁷

¹

-

3.3.4 数値計算結果

数種のβ に対する振動波形の計算結果を図3.8に示した これらの振動波形は ρ→1 (準線形微小振動〉の場合は正弦波へ， また0→o (強非線形 ， 超低周波 振動〉の場合は振幅1の矩形波に漸近していく. 図3.9は， 振I�b{ Axおよび謝波振 幅Bの周波数応答を示した. 本節に限り調波振l幅Bの定義は次式による.

X(T) =

'LBn sinnm (3.3.45)

n註l(odd)

ρ→1 では， 振幅および各調波振幅は，

Ax →2

-}

2(1-f2)/3→0，

ì

Bn→4J2

�

[(1-ρ)/12r→0， (f2→1 )

j

なる関係が存在し， 振幅， 振動 数ともに零 に漸近し， またAx →B}となる J主に

。→0では， 振幅， 調波振幅ともに非零の定数に漸近する. すなわち，

(3.3.46)

十ふ

の

そ ，刀

る

あ

で

h天 円楕 半四 hツ

お

ん」

の

門MM，年ハ

周

十Am ---J m「

ハU

X

3，→

Jl

do

--x

r(A

W久 生

L十仇情

→了

ん

→

澱

AL

B

振

幅 の 周

(3.3.4

7)

2 (β-f2，，1

X l ⁺

I ... �

^{:: 0}

^I

^{= 1 ， ただし，}f20 <ρ< 1， ^ì

人1-f20 ) ^{V '}

I

Ax

=

1，ただし， O<f2<f2o'

�

ここに， f20 < 0'""1/2

I

なる方程式でよく近似されることがわかる.

図3.10には，各種特性値と母数ノ守ラメータμ。 および補母数μ1との関係をぶし た μ。 →0においては，

(3.3.48)

Ax→

);;

_，

₁

1-f2→3μ。/4J 方， μ1→0においては，

(3.3.49)

l-Ax →μ1/4，

ì

0→π/ J21n(16/ μJ

J

なる漸近特性を有する.

(3.3.50)

ウU勺I

1

π 2π

Independent variable .!万

台、___"

N

.8 ロ 0

コ

∞ 。

1i

図3.8 対称ソフトスプリング系の解の波形

吋、〕勺/

1 1-.0

0.1 0.01

甘

5H

^{ω 口}

巳日 cu

4

1

国

に g 5

^J

_0.1

で同

8Z

^{ロEL I}

0 01

0.001

5 nu

崎区

1 1.5

π包^ロ

且5

コ 4

に'0."-

‘ _(/)

10 'U ^4A

、、 _... ヨ

、‘

国

に D

_同

同E

J

0.1

_‘ _.圃圃._-園町・・・・旬圃圃.Ax _{-、 句、}

1

旬、\

、、ミ

.、

、旬 ^圃‘

国

4/nπ←B_ (n ;odd) ^.、

5

²

0.01 0.5守

ミ?に

""0 ^ω

0.001卜|

I I I I I tI ^{‘" "" "}

，寸一一

。

z

0.01 0.1 1

^。

0.5 1

Frequency

.0

^Frequency

.0

図3.9 対称ソフトスプリング系の振幅の周波数応答

-74-

1 0.1

Complementary modulus parameter μ1

10-40 10-20 0.001 0.01 10-60

10

1

同町二貯dhqlHh

A布、で1HqhH/ / /

JO

...05'"

10 「

1 0.1

0.01

0.001

HFhsdhqlH

hq

hH布、lHRHで

0.001

v、/f \

J4U\

。\\

一

\ 噌i

\ \ \ 、

0.1 0.01

Modulus parameter

1 0.1

0.01

0.001

μ。

補母数料 と諸特性値の対応 母数μ。，

図3.10

戸、JづI

3.4 対称スナップスルースプリング系両振りモード

スナップスルースプリング系の運動方程式は次式で表される.

d 2X

-dT7-X+X" 3 二0 ， 周期(T)二2π/Q (3.4.1)

この 系には， まったく性質の異なる二つの独立な周期解(“両振りモード"および

“片振りモード"の解〉が存在し， どちらの モードの解が実現するかは初期条件に 依存して決定される. 式(3.4.1)で表される系には 1個の不安定平衡点(X = 0)と 2 個の安定平衡点(X = i=1)が存在し， 両振りモードの解は， 不安定平衡点である変 位の原点(X = 0)を中心として， 3個の平衡点をすべて通過する形で振動を行い，

その 振幅Axは， IAxl> J2である 一方IAxl<nとなる場合には， 片振りモード の解となる. 本節では， 両振りモードの解の取り扱いについて述 べ， 片振りモー ドについては次節にて論じることとする.

解は対称ノ\ードスプリング系と同様cn関数にて記述できるが，母:数の値は，

1/2 <μ。<1 の範囲に限定される. まず，解を次の ように仮定する.

X(T) = Axcnu(T)， Ax = max(T)IX(T)1 ;定数

式(3.1.18)および表 3.2 を参照することにより次式を得る.

AX2 二2μ。/(2μ。-1)，ただし， 1/2 <μ。<1

(3.4.2)

(3.4.3) (KOQ)2 ⁼ 1/ (2μ。-1)，ただし， 1/2 <μ。<1 (3.4.4) 母数の 変域が 1/2<μ。<1 であるので， 2.1 節の議論に従い， 式(3.4.3)， (3.4.4) を

μ1を用いて変形 すると，

1-^，u

Ax" = 2一一:l(=1-2μ1 2�∞付μ1二0�1/2) (3.4.5)

(

^一一一

叫 r

^{. 1 = =} _1-^�

_�

^{� (= 一}1 �∞付μ1二0�1/2) 2Qノ 1-2μl

両式から， 陽なμlを消去すれば次式を得る.

1+1 -竺左手五| 二AX2(=2�∞付μ1 = 0�1/2)

\. 2Qノ A

(3.4.6)

(3.4.7)

3.4.1 強非線形低周波振動(Ax→，)2， Q→0， μl→0)

1)振幅Ax(^< J3→，)2) が与えられる問題

-76-

式(3.4.5)を変形して， 与えられるAxに対して， μ1を得る表示は次式となる.

μ1 = (A x 2 /2-1) / (A x 2 -1)

= (A x / 12 -1)( 2 + A x / J2 -1) / (A x 2 -1) (3.4.8) 上式 からAxを与えて直接にμ1を計算しよう とするとき，

<<μ1/ AX)) = 6/ (2Ax 2 -1)(Ax2 -2 ) →∞， (AX→，)2) (3.4.9) であるので，Ax →12 においては，Axの誤差が異常に拡大されてμ1 の誤差にな る恐れがあるため，Axではなく (Ax-12)相当値を入力データとするよう注意す る必要がある. そのため， 式(3.4.5)を次式のように変更する.

C(2+ C) =μ1/ (1-2μJ，

ì

ただし，C三Ax/12-1

j

これから，

μ1二α/ (1 + 2a)， (α →0+)

ì

ただし，α 三C(2+C)， (C→0+)

j

上式によりμ1 μJα (C)Jを求める際の誤差増幅率は，

(3.4.10)

(3.4.11 )

<<μ1/ C)) = <<μ1/α ))<<α/ C)) = [1/ (1 + 2a)][2/ (2+ C)]→0(1)， _ì

ト (3.4.12)

(α→O，C→0， Ax →12)

J

であり， Ax →J2において安全となる. μ1 からqlを得るには， 式(2.1 .15)を試 行解法により解くか， ま たは 2.4節で作成したql =qj(μJ なる|場表示多項式[豆、

(2.4.6) ]を用 い る . さ らに ， 式 (2.1 .13) か らQ=Q(qJを， 式(2.1 .14) かり

�=民(qJを得る. これら一連の計算に関しては，2.5節より誤差の増幅に関し ては安全であることが確認済みである.決定されたμlヲQ，�を式(3.4.6) に代入す れば0が特定される. ここで， この0の計算においては， 式(3.4.6) ではなくず (3.4.7)を利用すると， 与えられたAxが部分的に関与するから若干の精度向上か 期待できるだけでなく， 形も簡潔であるから， より好ましい表示といえる.

-77-

→

(

Ax / K1(μ1 =1/2) (→∞) π/川x/vÍ2-1)]

(�O)

^;Ax→;Ax →∞ l

� r

(3.4.13)

上式の誤差 増幅率を計算しておくと，

((f2/ Ax)) = Ax2 /(Ax2 -1)

+((Q /μ1))((μ1/ AX))+((ζ/μ1))((μ1/ Ax)) (3.4.14) であり，μ1 二，μl(Ax)を計算する際， 式(3.4.11)を用いる よう注意さえすれば，

((f2/Ax))→0(2)， (Ax→vÍ2) (3.4.15)

であり，誤差の増幅に関しては問題ない. こ れらf2，Qによって式(3.1.15)により フーリエ係数が決定される.

2)振動数。(< 1.03→0)が与えられる場合

式(3.4.6)を用いて， 与えられたρからμ1 あるいはq1を計算するとき の誤差 の増幅率を計算すると， 次の結果を得る.

((θ/ q)) = ((Q / q)) -((1(， / q)) + ((μ1/q1))[μ1/ (1-2μJ]，

ì

→O(((Q / q)))二0(11/1nμ1)1 →

�

^{� ( ， 5μ1→0)}

{

l

^{∞， (μl→1/2)}

J

^(3.4.16)

ゆえに， 111→0では与えられたβ からq1 あるいはμ1 を得るときに危険であ り，逆にμ1→1/2では，μ1 から0を得るときに危険となる.そこで，対祢ソフ

トスプリング系の場合と同様に， 式(3.4.6)を変形 して次なるL1F を定義する.

L1F二1rQ / 2Q -1 = 1/ K

f1三五 -1

^{= fun.(qJ， ì}

.・.πβ/2Q= 1+ L1F，

ト

π2 /2Q = (π/Q)(l+L1F) J

他方， 式(2.1.13)から，

q1二exp(-π2/2Q)=qoqndF 1 ただし，qn三exp(-π/f2)，

ト

qρðF三exp(-^7[ L1 ^F / f2) )

ここで，L1F = 1/ 1(，f1弓瓦-1 の情報落ちを防ぐ形に変形 すると，

(3.4.17)

(3.4.18)

A ニ8�ql(1+bJ4 - � (1 + 2�)

l

K]

$三五

^(1+KJ

$三五 ↓

= fun.[qJ( L1F， Q)] = O(12qJ， ql→o )

-78

(3.4.19)

与えられたβ に対して， この式をL1Fに関して試行的に解くことができる. 式 (3.4.19)の右辺のql中の4を仮定し計算すると， それに続く改善されたL1F を 得るからL1Fの試行解の出発値はL1F二0 とおけばよい L1F = L1F( Q)が得られ ると式(3.4.17)よりQ= Q( Q)が求まって， 式(3.1.15) よりフーリエ係数も定ま る. このQ=Q(β)の計算に関しては，ここで必要とされる変域[μl二0�1/2 ，

(.(2= O� ∞) ]のすべてにおいて， 精度的に安全であるこ とが確認済みである.

ここで，μ1-μl(ql)の計算の際必要となるqlを特定しようとするときには若工 の注意が必要となる. 図2. 7より，ql→0 ではqlから μ1 を計算することは精度 的にはまったく問題ないが，。あるいはQからqlを計算する際には対称ソトス プリング系の場合と同様の問題が生ずる. すなわち，οから式(3.4.18) によるに せよ，あるいは一旦求められたQから式(2.1.13) によるにせよ，。→0 において ql を求める場合に， その誤差 増幅率 を計算してみ ると，

((ql / / .(2))→π/.(2→O(11nOI)→∞， [式(3.4.18)J _ì ((ql / /Q))→ポ/2Q→π/.(2→O(11nOI)→∞， [ 式(2.1.13)J，

ト

ただし，ql→o ^j

(3.4.20)

となり計算精度上問題が生じる. 実 際問題として，計算機内部の最大~最小値は，

1075 �10-75程度であるから，qlの最小値 を10-75とすれば， それに見合うοある いはQの最小値ならびに式(3.4.20)の誤差の増幅率は，

。> 0. 02， Q > 0. 03， ((ql / /θ)) < 157

となって， 単精度 計算以外ならば， 何とか許容できる値と考えられる. しかし，

.(2，Qがこの値以下になると， 計算機内では事実上。=0と同等に扱われてしま うこ とに注意する必要がある 以上より得られた μ1=μJ(qJ を式(3.4.10) に適用 すればCが求められる. すなわち，

C二 μ1/ (1-2μJ

1 +

�

(1- μJ/ (1 -2μJ (3.4.21 )