ショート・ペーパ

計測自動制御学会産業論文集
Vol.9, No.22, 163/166 (2010)
163

片持ち梁の固有振動数計測に基づくき裂位置の同定
中住昭吾*，西郷宗玄**
Crack position identification based on measuring eigenfrequencies of cantilever beam
Shogo Nakasumi*, Muneharu Saigo**
The crack position in a cantilever beam is identified. The vibration experiment is conducted and the eigenfrequency data
is obtained by measuring the acceleration of the beam tip. The extended finite element method is used to analyze thevibration problem. The evaluation function is calculated using the eigenvalues of both experimental and numericalsimulation, and the identification is conducted based on the gradient of the evaluation function. The accuracy of longitudinal
direction is much more higher than that of thickness direction. That difference is caused by the characteristic distribution ofeigenvalues.
Key Words. Identification, Crack, Eigenfrequency, Nondestructive Inspection

y１. はじめに構造物のき裂・空孔等の損傷箇所を同定する非破壊検査とx ,y )x (mmCrack (Real)

Experiment& measuring

λim(a) Actual crack position in actual beam member(x , y )Crack (Estimated)

して超音波やマイクロ波を用いた検査法が用いられるが，こy

れらは構造物表面をスキャニングするための作業量コストの(b) Estimated crack position in calculating scheme

exe

Solving eigenvalue problem

λie多さが課題である．一方．構造物の動特性，すなわち固有振Fig.1 Actual and estimated crack position

動数に着目した同定方法が近年注目されている 1)．この方法

裂端座標はそれぞれ( x m , y m ), ( x e , y e )と表記される．この概はスキャンのための労力を大幅に軽減できることが期待でき念図を Fig.1 に示す．同定の指標としては，実在き裂と推定る．そこで本論文では，片持ち梁の固有振動数計測実験を行

き裂の固有振動数の差を入力値とし，両者が一致するときにい，その計測結果を基に数値計算によるき裂同定を行なう．

最小となるような値を出力する関数（以後，評価関数と呼その際，拡張型有限要素法と呼ばれる手法を用いることにすぶ）を設定する．評価関数式は文献 1)と同様に(1)式を用いる．

る．

(1)式は振動モードにより複数個存在する固有値を平等に扱

2.固有振動数によるき裂同定の定式化とアルゴリズム

い，それら全体の誤差を最小二乗法的に算出するものである．

き裂同定の構造部材として片持ち梁を対象とし，き裂は表同定に用いるモードの総数である．本論文では n=3 とした．

ここで，固有値 λ i の右下添え字 i はモード次数であり，n は面から垂直に入るものとする．これよりき裂同定は長手方向位置とき裂深さの 2 個のパラメータで一意に決定出来る．実在の梁にその位置が未知であるき裂が存在し (計測により固有振動数は取得できる)，固有振動数情報からなる何らかのJ={(1 n∑ λ ie − λ im
2 i =1)λ}mi

2（1）

逐次変化する推定き裂位置に対応する λ ei を得るためには，指標に基づいて，計算上で仮定したき裂位置を移動させき裂

その度に有限要素法による固有値解析を行なう必要がある．

を同定するというプロセスをとる．以下では便宜上，前者をそれにはき裂位置変化を反映した有限要素メッシュ（以降メ実在き裂，後者を推定き裂と呼び，両者に関する物理量を右ッシュと呼ぶ）を用意する必要がある．しかし一般にメッシ

肩の添え字 m ならび e をつけて区別する．これに倣うとき

ュの作成は解析者による手作業の修正を必要とし，これが処* 産業技術総合研究所先進製造プロセス研究部門
茨城県つくば市並木 1－2－1 つくば東** 東洋大学理工学部 埼玉県川越市鯨井 2100
* Advanced Manufacturing Research Institute, Advanced
Industrial Science and Technology (AIST), Tsukuba-shi, Ibaraki** Faculty of Science and Engineering, Toyo University,
Kawagoe-shi, Saitama(Received May 18, 2010)

理の自動化の障害となる．本論文ではこの問題に拡張型有限要素法 (eXtended Finite Element Method : XFEM, 以下 XFEMと呼ぶ) と呼ばれる計算手法 2) を用いて対処する．このXFEM では，要素内での近似関数分布を一次多項式等の単純な関数のみに限定せず，より一般化された形を許容する．こTRIA 0022/10/0922-163©2010 SICE

164

SICE T. on Industrial Application Vol.9 No.22 2010

JJJJHeaviside functionC

C Enriched node with

CC

L=500mm
H=6mm
B=6mm

J Enriched node with

Crackr yθ
CYX

Acceleration sensor(B&K4374)

SupportZ

れを近似関数のエンリッチ (Enrich)と呼び，これにより要素内部にき裂をはじめとする複雑な物体形状をメッシュ形状に依存せず構築できる．以下では，2 次元モデルに対してメッシュとき裂が Fig.2 に示すような位置関係にある場合の
XFEM によるき裂周辺の変位場の定式化を述べる．き裂周辺に位置する節点に対して，それに結合する要素がき裂によって完全に切断される場合には J 属性（変位の不連続性を表わす），あるいは要素が完全には切断されず要素内部にき裂端を含む場合には C 属性（き裂端の応力集中を表わす）と表記する属性が生じるものとし，それらの節点には上述したエンリッチの関数に対応する節点自由度が新たに付加される．
その結果，き裂周辺に位置する要素 e 内の変位場 u e は(2)式で表わされる．

∑L

E=68GPa
ρ=2700Kg/m3Fig.2 FEM mesh and configuration of enriched node

ByAluminium beam

4ca

singular function
Normal nodexV1a=3mmc=1mmzHxFig.3 Overall view of experimental system
振動解析では(2)式より全体剛性行列 K と全体質量行列 Mを導出し，(5)式の一般固有値問題を解き固有値 λie を得る．

&&+ Ku = 0Mu(5)

固有値 λie が求まるごとに推定き裂位置を更新するが，その際メッシュの修正は不要であり，(3)式, (4)式のエンリッチ関数の更新だけを行なえばよい．
(1)式の評価関数値 J の最小化には最急降下法 4)を用いる．
すなわち x 方向移動については(6)式となる．ここで左上添え字 i は計算ステップ数であり，また α は１ステップ当たりの移動量に関する比例定数である．(6)式の ∂J ∂x は(1)式より(7)式となる．(7)式の ∂λie ∂x は固有値 λie の感度算出式で

∑∑4

∑

ue =
NI dI +NI H ⋅ bI +N I ψ K ⋅ c IK
I =1I ∈JI ∈C
K =1(2)

ここで右下の添え字 I（I=1～4）は要素内で定義された節点番号であり N I はその節点に対応する形状関数である．通常の有限要素法では(2)式を右辺第一項で打ち切るが，J 属性または C 属性の節点についてはそれぞれ第二項または第三項が加わる．また， H と ψ K ( K = 1 : 4 ) は上述したエンリッチ関数であり，それぞれ(3)式，(4)式で表わされる．前者はあり(8)式で求まる 5)．(8)式の u i は固有値 λie に対応する固有ベクトルであり，感度行列である ∂K ∂x と ∂M ∂x は差分近似により導出する．なお(6)式~(8)式は y 方向についても同様の手続きとなる．
i +1x = i x − α ∂J ∂xn(6)(λ − λ ) (λ )∂λ ∂x = u ( ∂K ∂x − λ ∂M ∂x ) u∂J ∂x = ∑  ∂λie ∂x Teimimi2(7)i =1eiTieii(8)ヘビサイドのステップ関数，後者は線形破壊力学により導かれるき裂端変位場の基底関数である．(3)式 y と(4)式の r , θ３．振動数計測とき裂同定

はき裂端を原点とするそれぞれ直交座標系と極座標系の座標片持ち梁の振動数計測実験における部材寸法，物性値を値である．(2)式右辺の d I ， b I ， c KI はそれぞれ対応する節

Fig.3 に示す．部材はアルミニウム合金である．梁長さは点自由度ベクトルである．

50cm で一端を万力で固定し，梁先端には微小一軸加速度計

H =1(y ≥ 0), −1 (y<0)(3)

xy 面内の２次元梁曲げ振動を計測することになる．き裂をψ K = [ψ 1 ψ 2 ψ 3 ψ 4 ] θ= r sin
 2cosθ
2sinを取り付けた．計測される加速度は y 軸方向のみとなるため，幅 1mm，深さ 3mm の矩形ノッチで表現しその位置が固定端θ2cos θθ(4)cos cos θ 2

より 100mm または 350mm となる場合を設定した．上記のき裂を有する場合，およびき裂がない場合の計３ケースにつ以上が拡張型有限要素法における変位場の表現方法である．
また，き裂位置により定まるエンリッチ節点の選定アルゴリ

いて，梁をハンマで殴打し伝達関数から固有振動を計測した．
数値計算に用いるメッシュの分割数に対して固有振動数が）一十分に収束することを確認する必要がある．各要素を正方形意には定まらない．より詳細な定式化や理論的背景は文献とし，厚さ方向の要素分割数を 1～6 とした場合に算出され2)3)る固有振動数を Table1 に示す（き裂がない状態）．ただしズムは種々の幾何学的計算方法があり（たとえば文献等を参照されたい．

3)

Normalized
frequency

Mode1
Mode2
Mode3

Number of elements in thickness direction
123
4
56
1.219
1.059
1.027
1.015
1.010
1.007
1.219
1.058
1.026
1.015
1.009
1.006
1.218
1.057
1.025
1.013
1.008
1.005

Table 2 Eigenfrequencies of cantileverbeam
Experimentfrequency[Hz]

No crack 100mm
mode1
mode2
mode3

18.4
116.4
325.5

17.1
113.6
313.2

350mm
18.2
110.8
304.0

Numerical
No crack No crack 100mm
(Y68)
(Y60)
(Y60)
19.7
18.5
19.4
123.2
115.7
126.3344.6323.7
349.3

350mm
(Y60)
20.1
123.8
339.7

ベルヌーイ・オイラー式の理論値にて正規化した値である．
この結果，厚さ方向に 4 要素（長手方向には 333 要素）でほぼ十分収束したと見なし，以降の数値計算でこのメッシュを用いた．また上述した実験条件より，片持ち梁の変形が大き
い固定端付近ではｚ軸方向の変形拘束が作用すると考えられることから，(5)式の K は平面ひずみ状態とした．実験計測
および数値計算より得られた振動数を Table2 に示す．

6Real cracktip

4200x=109mm
y=199mm
J=7.49x10-4

ち消すために数値計算で使用する材料の Young 率を微調整

200300Longitudinal location : x [mm]400500642S.P.1

S.P.3

Real crack

Converged pointx=181mm
y=193mm
J=1.58x10-3Real crack tip

S.P.200

Converged pointx=335mm
y=218mm
J=6.96x10-4

S.P.4

100200300
Longitudinal location : x [mm]400500(b) Crack position : 350mm from clamped end
Fig.4 Crack identification path
Eigenvalue (mode2)

機と計算との間に生じたモデル化誤差と捉え，この誤差を打S.P.4

S.P.2

100

Eigenvalue (mode1)

材料の特性となったことが原因と考えられる．この誤差を実Converged pointx=315mm
y=291mm
J=1.39x10-2

Converged point

ラムにて解析した．その結果計算値の方が 5.9%～7.1%程度
境界条件が実験では厳密に実現されておらず，より軟らかい

S.P.3

S.P.1

Real crack

まずき裂がない状態で固有振動を本論文で使用するプログ計測値よりも大きい値となった．これは主に片持ち梁の固定165(a) Crack position : 100mm from clamped end
Thickness location : y [mm]Table 1 Normalized eigenfrequencies

Thickness location : y [mm]計測自動制御学会産業論文集 第 9 巻 22 号 2010 年

422001004300 400

2100300 400

Eigenvalue (mode3)することとする．Young 率を 68GPa から 1GPa ずつ下げて固200

有値解析を繰り返した結果，実機材料よりも 12%程度小さい値である 60GPa で計測値とほぼ同等の振動数となることがわかった．よってき裂がある場合に対応する数値計算には
この Young 率 60GPa を用いた．Table1 の Y68, Y60 の記号は42

これらの違いを意味する．実機と数値計算間のモデル化誤差

100はこの Young 率の修正により第一次近似的には補正された

200300 400

Fig.5 Eigenvalue distributionと考えられる．しかしこのような補正を経た後でも，き裂が存在する場合には計算値の方が計測値よりも振動数が全体的ような場合でも収束地点の評価関数値を直接比較することでに大きい値となり，実機の方がより軟らかい特性を示すことどれが大域的な最小値であるか判断することができる．

がわかった．これが同定の精度に及ぼす影響は後程考察する．

つ ぎ に 同 定 の 精 度 を 考 察 す る ． き 裂 位 置 (100mm と実験計測した振動数を基に(6)式~(8)式のアルゴリズムにて350mm) の違いによる影響は見られなかった．また，実在き同定した結果を Fig.4 (a), (b)に示す．試行錯誤の結果，(6)式

裂端と収束した推定き裂端の座標値のずれ量 (x 方向，y 方の α は x,y 方向にそれぞれ 1× 10 4 ， 1× 10−4 とした．Fig.4 は向) をそれぞれの方向の梁部材寸法値で除した値は，Fig.4

推定き裂のき裂端座標の移動履歴 (同定経路履歴) を示して(a), (b)について長手方向が各々1.8%, 3.0%, 肉厚方向がおのおおり，座標値がほぼ十分収束したと見なせる時点で計算を終の 16.8%, 13.7%となった．前者については非常に高精度な同了した．おのおのの収束地点での座標値 (x,y) および(1)式の定が行なえた．一方，後者については前者ほどの精度は得ら

評価関数値 (J) を図中に示している．Fig.4 は(a), (b) 共に 4 かれなかった．

所ずつ初期値 (Start point) を設定しており，図中にそれぞれ同定誤差に方向別の大きな相違が生じた原因を考察するた

S.P.1~S.P.4 の記号で位置を示している．同定プロセスを経て，め，評価関数の入力値である固有値の変化を調べる．き裂の経路は最終的に 2 地点に収束した．これは評価関数が複数個存在により固有値は一般に減少するが，その減少度合とき裂所の局所的停留値を持つ多峰性分布となることによる．この位置の関係を調べるため，き裂位置を梁全体領域で変化させたときに数値計算で求まる固有値( λ1 ~ λ3 ) の分布を Fig.5 に166

示す．長手方向変化についてはモードごとの固有の減少パタ

長手方向位置が非常に高精度に同定できた．今後はき裂深さーンが明瞭に読み取れる．一方で深さ方向変化についてはどの同定の精度を向上させることが課題である．のモードもき裂が深くなるにつれ単調に固有値が減少する．
これらの理由により，計測値が計算値よりも全体的に軟らかい特性を示しても，長手方向位置は複数の固有パターンから最小二乗法の原理によりある程度正確な位置を同定でき，一方深さ方向にはどのモードも同じ傾向を示すため，正確な値を検出することが出来なかったと思われる．

４．結言
実験計測から得た片持ち梁の固有振動数情報を利用し，数値解析上のき裂位置を連続的に移動させて振動数誤差最小原理に基づくき裂同定システムを構築した．同定の精度を検討し，肉厚方向はややき裂が深めになる傾向となったものの，

参考文献
1） 堀辺忠志，浅野直輝，岡村弘之：き裂を有する弾性ティ
モシェンコはりの固有振動および GA に基づくき裂同定，
日本機械学会論文集(C 編), 66, 1756/1761(2000)
2） N.Moes, J.Dolbow and T.Belytschko : A finite element method
for crack growth without remeshing, International Journal for
Numerical Methods in Engineering, 46, 131/150(1999)
3） M.Duflot : A study of the representation of cracks with levelsets,

InternationalJournalfor Numerical

Methods inEngineering, 70, 1261/1302(2007)4） 田村明久，村松正和：最適化法，p.102，共立出版 (2002)
5） 山 川 宏 ： 最 適 化 デ ザ イ ン ， p.82, 培 風 館 (1993)