68-１
フーリエ変換による緩和スペクトルの解析手法の構築
―木材の粘弾性評価に向けた試案―
明瀬 航 1. はじめに 木材の粘弾性は, JIS 規格や構造計算の中で考慮さ れる程, 広く知られる特徴的な性質でありながら, そ の評価方法を確立する上では多くの課題が残る. 木材 に固有の問題も複数あるが, その理論的枠組を与える レオロジー分野でも共通する問題として, 評価尺度の 確立が挙げられる. 一般に粘弾性特性は, 対象となる 物体の性質を線形モデルの重ね合わせを軸とする, 緩 和スペクトルや遅延スペクトルに置き換えて解釈する. これらのスペクトルは, 一般に動的試験から求められ, 温度－振動数等価性を利用することで1), より広い範 囲の性質を把握する試みが成されている. この方法は, 試験設備があれば有用であるが, 木材の試験を前提と する場合, 温度や水分量との関係が明らかでない他, 寸法効果の影響も未解明なため, 試験が実施可能な設 備は我が国にも多くない. これらを背景に, 本研究で は静的試験より緩和スペクトルを求める方法の確立を 目的に, 数値計算上の課題の解決方法を提案し, その 妥当性を数値実験により検証した結果を述べる. 2. 緩和スペクトルの計算方法 Renardy が示したフーリエ変換と緩和スペクトル の関係性を参考に2), まず連続関数系での計算方法を 整理する. 緩和スペクトルは次式で定義される. (1) 式 (1) で, G : 緩和弾性率, F : 緩和スペクトル, t : 観測時間, τ : 緩和時間である. ここで, 緩和時間 τ を 広範囲にとる際に便利が良いため, 対数として整理し 直すと, の関係から (2) となり, 式 (2) で H : 緩和スペクトル と定義され, ここで緩和スペクトルと緩和弾性率の関係が導かれる. 続いて, とおき, 変数変換を行い, さら に, 各関数を (3-1, 2) と置き換える. 式 (2) に式 (3)を代入する. H(lnτ) の 値域が正, 積分変数を であることを踏まえ, 積 分範囲を一般化してある. Φ をフーリエ変換すると, (4) と表せる. 式 (4) に x = u - σ を導入する. (5) 式 (5)で s = ex を導入し, 積分変数を置換することで, (6) が得られ, Γ 関数を導入して表記すれば, (7) となり, Γ 関数を介したフーリエ変換後の緩和弾性率 と緩和スペクトルの関係式を導くことができた. 以上の関係式を, 実際の解析手順で計算式を並べた フローチャートを図 1 に示す. 図 1 右に緩和弾性率 から緩和スペクトルを求める本論文の目的とする手順 (以下, 正手順), 左には緩和スペクトルから緩和弾性 率を求める手順 (以下, 逆手順) を示す. 図 1 解析のフローチャート 式 (2) より, 緩和スペクトルを数値積分することで, 容易に緩和弾性率を求めることができ, 逆手順は示す 必要はないが,フーリエ変換による緩和スペクトルの 計算方法は前例がないため, 周波数領域の解析値の妥 当性を検証することができない. そこで, 梗概では紙 面の都合上省略したが, 事前に逆手順の解析の特徴を 把握した上で, 次章で行う正手順の解析を行っている. σ Φ Φ Relaxation Spectrum→Relaxation modulus of elasticity ( Inverse procedure )
σ σ σ
τ τ τ
Relaxation modulus of elasticity
→ Relaxation Spectrum ( Regular procedure ) σ Φ Φ + Auxiliary Section Φ Γ Φ DFT IDFT Φ Φ Γ σΦ DFT IDFT Noise Processing Continuity Revision
68-2 3. フーリエ変換による緩和スペクトルの解析 本章では, 実験データから緩和スペクトルを解析す ることを目的とするため, 離散化の影響が付きまとう. 実験で得られる緩和弾性率の離散データが時間幅 TS, 時間間隔 Δu でデータ数は N = TS
/
Δ + 1 となる. この離散データのフーリエ変換後の周波数領域にお いて, ナイキスト周波数 fN/2, 振動数間隔 Δf は式 (8) の通りに定められる. (8-1, 2) 以上の離散式に基づき, 本章の解析を行う. 3.1 周波数領域における誤差とその解決法 前章に示した計算手順で解析を行うと, 緩和スペク トルの解析結果は非常に煩雑な結果が得られてしまう. その最大の原因として, 周波数領域における誤差にあ る. その誤差を引き起こす主な要因は 2 つある. 1 つはフーリエ変換による周波数領域におけるノイ ズである. ここでノイズの存在を確認するため, 式 (9) を用い, 周波数領域の理論値と解析値を比較する. 𝑓(𝑥) = 𝑒− 𝑥2 2𝑎2⟺ 𝑓̂(𝑣) = √2𝜋 𝑎 𝑒−𝑣 2𝑎2 2 (9) 図 2 (a) にはフーリエ変換前のガウス曲線を示す. 図 2 (b) にはフーリエ変換後の理論値と解析値を示す. 解析条件として, TS : 50, fS : 1000, 標準偏差 a : 1.0 と 設定し, 解析を行った. 理論値と解析値はほぼ合致し ているように見えるが, 拡大図に注目すると, 0 付近に 微小なノイズが発生していることが分かる. 図 2 フーリエ変換前後の理論値と解析値 ノイズの原因として, 離散データに対しての演算, 窓関数のサイドローブの影響等と推察している. 解析 プログラムでは, 倍精度実数を使用しており, 時間幅, 時間刻みを変え, データ数を増加させてもこの微小な ノイズが消えないことも確認している. 次に, 図 3 (a) に, 図 2 (b) のフーリエ変換後の関 数, C : 0.1 の Γ 関数の実部, Γ 関数除算後の関数, それ ぞれの理論値と解析値,計 5 つの関数の絶対値を縦軸 に対数軸をとり, 示している. 縦軸に対数をとってい るため, グラフが下に行くほど値が 0 に近づいてい る. Γ 関数除算前後の関数の理論値と Γ 関数は周波数 が増えていくにつれて, 0 に漸近していく様子が分か る. 図 3 ガンマ関数除算前後の解析値と理論値の比較 ここで, 理論値と解析値を比較することで, フーリ エ変換によって v : 8.80 から明瞭にノイズが始まって いることが分かる. フーリエ変換によるノイズは振幅 が比較的一定となるため, 対数軸をとると理論値と外 れ横に伸びている状況が見て取れる. そのノイズが現 れた関数をΓ 関数で除算することで, ノイズが始まっ た点から値が徐々に増えていく. 図 3 (b) には Γ 関数 除算後の関数の理論値と解析値を示す. 低周波数の領 域では理論値と解析値が一致している様子が分かる. 理論値では v が 10 を過ぎたあたりから 0 に収束し落 ち着いていく. しかし, v が 20 を越えたあたりから解 析値のグラフが縦軸の正負に絶対値が膨大に増えてい る. そこで, の全データの最小値となる点からノイ ズが始まると考え, その時の周波数を遮断周波数 vh とし, vh 以上の周波数の のデータを全て0 にし, 演 算を実行することで誤差を最小限に抑える. もう1 つの周波数領域においての誤差を誘発する要 因が, 関数の不連続性である. 有限フーリエ変換は, 無限に繰り返し循環している波の一周期分を取り出し, 計算しているに過ぎない. つまり, 関数の始点と終点 が一致していなければ, フーリエ変換後の周波数領域 で好ましくない影響があることは容易に想像できる. 本計算では, 周波数分解能が重要であると考え, 窓 関数は矩形窓を採用したが, 関数の不連続性を放置し てしまうと, 正確な緩和スペクトルを計算することが 出来ないと判断した. そこで, 連続性を確保するため に, 例として, u : -25 ~ 25 の時間幅で解析を行うなら ば, 式 (10) の演算を利用することで, 平行移動量 p を求め, 各時刻の緩和弾性率に足すことで, 平行移動 を行い, フーリエ変換前の関数である Φ の連続性を保 つことができる. (10) 図 4 緩和弾性率の平行移動と連続性確保の様子 0 1 -5 -3 -1 1 3 5 f (x) -1 0 1 2 3 0 5 10 15 10-15 -10-15 8 10 12 (a) (b): Theoretical Values : Analytic Values
0 5 10 15 104 1 10-4 10-8 10-12 10-16 10-20 vh = 8.80 min -100 -50 0 50 100 0 5 10 15 20 25 30 vh = 8.80 (a) (b)
: Theoretical Values : Analytic Values
-25 -15 -5 5 15 25 G u -25 -15 -5 5 15 25 u p p e-25C{G(-25)+p} e25C{G(25)+p} (a) (b)
68-3 3.2 仮想の時間幅に対する緩和スペクトルの解析 本節では, 前節で提案した解析手法がより複雑な緩 和スペクトルに適応できるか検証した結果を示す. モデルの概要は, 表 1 に示す. 式 (11) を用いるこ とで, 図 5 下段に示す case 1 ~ 4 の緩和スペクトルと なる. 緩和ガウス時間 Iσ とは正規分布の平均であり, 緩和ガウス頻度 IH として係数を与えることでより自 由な曲線を描けるように設定した. 表 1 緩和スペクトル曲線のモデル概要 (11) case 1 ~ 4 の緩和スペクトルから緩和弾性率を作成 し, フーリエ変換による計算結果と事前に解析した逆 手順の解析結果との比較を行う. 解析条件は TS :50 , fS : 100, C: 0.1 とし, 解析を進める. 図5 に case 1 ~ 4 の解析途中を含めた結果を並べ て示している. 上段に緩和弾性率を, 中段には の虚 部を, 下段には緩和スペクトルを示している. 各 case でフーリエ変換によるノイズが始まる点までは逆手順 で解析した値を辿っていることが分かる. ノイズが始まる点と遮断周波数 vhが外れているの はcase 1, 3 であり, 解析結果にノイズが生じ, 負側端 部で拡張している様子から, ノイズ処理が解析結果に 大きく影響を与えることが分かる. case 2, 4 は遮断周波数 vhとノイズが始まる点が一 致しており, case 2 の解析結果はモデルと良好に一致 していることが分かる. しかし, case 4 では緩和スペ クトル負側端部にノイズではない縦軸に引き伸ばされ ているような曲線がある. これは逆フーリエ変換後の 関数に対して, リンク効果 3)が影響したものと考えて いる. 前述したように, 有限フーリエ変換は, 関数の 周期性に基づき解析しているため, 実関数が端部で不 連続であれば, 周期性を保つため, 結果をゆがめてし まう. リンク効果についての対策は次節に述べる. 3.3 応力緩和実験を想定した緩和スペクトルの解析 緩和弾性率のような緩やかに変化していく曲線に 対して時間幅が短くなることで式 (8)より周波数領域 における解析の振動数間隔が短くなることは, 解析結 果に悪影響を与えると予測できる. また, 前節で問題 提起したリンク効果によって緩和スペクトルの解析結 果をゆがめることが分かった. 得られた時間幅内で計算を行うことで引き起こさ れるリンク効果, 時間幅が小さくなることで振動数間 隔が短くなる, これら 2 つの問題を解決するために, 得られたデータの時間軸正負に補助区間として付け加 え, 解析を行うこととした. 図 6 補助区間設定とリンク効果の様子 図 6 に補助区間を設け, 解析を行う模式図を示す. 左上に実験で得られた緩和弾性率がある. ここで中央 上の斜線で囲まれた時間幅の緩和スペクトルを求めた いが, この時間幅では解析に不十分である可能性があ るため, 緩和ガウス時間 Iσ : -30, 30, 緩和ガウス頻度
case 0 case 1 case2
緩和ガウス時間 Iσ 0 5 -5 -10 0 10 -15 -7.5 0 7.5 15 緩和ガウス頻度 IH 1 5 5 10 3 10 3 5 7 10 2 標準偏差 a 5 1 1 要素数 J 1 1 1 case3 case4 1 3 3 5 -5 0 5 10 15 20 u G knnown Relaxation modulus of elasticity
-50 -30 -10 10 30 50
unknown Relaxation spectrum
Auxiliary section Auxiliary section
-50 -30 -10 10 30 50 u G
Relaxation modulus of elasticity + Auxiliary section -5 5 15 -50 -30 -10 10 30 50 u -50 -30 -10 10 30 50 + = link link Continuity Revision 図 5 case 1 ～ 4 の緩和スペクトルの解析結果 -1 5 -25 -15 -5 5 15 25 case 1 0 5 10 104 1 10-4 10-8 10-12 10-16 10-20 108 vh= 4.52 case 1 -1 5 -25 -15 -5 5 15 25 case 2 0 5 10 104 1 10-4 10-8 10-12 10-16 10-20 108 vh= 3.14 case 2 -1 10 -25 -15 -5 5 15 25 case 3 0 5 10 104 1 10-4 10-8 10-12 10-16 10-20 108 vh= 4.52 case 3 -1 5 -25 -15 -5 5 15 25 case 4 0 5 10 104 1 10-4 10-8 10-12 10-16 10-20 108 vh= 1.13 case 4 0 10 -25 -15 -5 5 15 25 G u case 1 0 50 -25 -15 -5 5 15 25 G u case 3 0 100 200 -25 -15 -5 5 15 25 G u case 4 0 5 10 15 -25 -15 -5 5 15 25 G u case 2
68-4 IH : 3, を時間軸正負に補助区間として設定し, 実験で 得られた緩和弾性率に滑らかな曲線を足すことで. 解 析条件の TS を 100 にし, 周波数間隔を短くする. か つ解析後の緩和スペクトルのリンク先を得られる緩和 スペクトルから外すことで, 解析結果へのリンク効果 による影響をなくすことができると考えた. 本節の解析で使用する時間幅は 20 で, u : -4 ~ 16 と設定した. これは 100 Hz の時間間隔で約 3 か月計 測した実験データの解析を想定した. ここで, 時間軸 負側は e-4 = 1.83×10-2 [sec], 正側は e16 = 8.89×106 [sec]となり約 100 日のデータとなる. 以下で, 実験を想定した時間幅の緩和弾性率のデ ータに対し,補助区間有無の解析結果を比較する. 対 象とする時間幅 u : - 4 ~ 16 で case 1 ~ 4 の緩和スペク トルモデルを切り取ったときに, 両端で顕著に不連続 性が生じる case 2, 4 を対象に解析を進める. 解析条 件は, C: 0.1, 補助区間無しで, TS :20 , fS : 100 補助区 間有りで, TS :100 , fS : 20 と設定し, データ数はどち らも2001 点としている. 図 7 に補助区間有無で比較した解析結果を示す. 左にcase 2 右の case 4 の結果をまとめて載せてある. 上段に の虚部を示す. この解析の理論値は分から ないが, 補助区間無しの遮断周波数はとても小さく, さらに振動数間隔はこれまでの解析条件よりも短いた め, 有効とみなせるデータは補助区間有りと比べてと ても少ない. 下段には, ψ を示しており, 補助区間無 しでは実験を想定した時間幅 u : -4, 16 で解析値がリ ンクし, 端部で値がゆがんでいるが, 補助区間有りで は, リンク効果が u : -50, 50 で影響するため, リンク 効果による解析結果のゆがみはない. しかし, 緩和ス ペクトルの解析値が補助区間と実験を想定した時間幅 の境目で急激な変化となるため, 解析値の端部に多少 誤差があるが, 緩和スペクトルの傾向を見るには十分 な再現度であると言える. 図 7 補助区間の有無による解析結果の比較 4. 従来の解析手法との比較 静的試験より緩和スペクトルを求める手法は, 以前 より検討されており, 例として, 木材の緩和スペクト ル曲線の解析を試みている桑村の論文 4)において使用 されているAlfrey- Doty の第 1 近似法 5)がある. 第1 近似法は, 緩和弾性率を微分することで, 緩和 スペクトルの近似値を得る. 図 8 (a) に case 4 の緩 和スペクトルの理論値とフーリエ変換による解析値, 第1 近似法による解析値を示す. 解析対象となるデー タが, 平滑な場合であっても, 第 1 近似法はピーク値 や起伏にわずかな差がある. また, 実験データでは, 本論文で使用したモデルで 得られる緩和弾性率とは異なり, 試験機の分解能や温 湿度の変化等のノイズが必ず含まれる. そこで, 振動 数0.01 ~ 50 の範囲で振幅 10 -4, 位相に乱数を発生 させ, 白色雑音を作成し, 解析対象の緩和弾性率に足 し合わせて, そのデータを用い, 解析を行った. 図 8 (b) に雑音を含んだ緩和弾性率から緩和スペ クトルを求めた解析結果を示す. 第 1 近似法は緩和弾 性率の傾きを求めているので, 雑音が解析結果に多大 な影響を与えていることが分かる. フーリエ変換によ る解析は平滑なデータからの解析ほど再現はできてい ないが, ノイズの影響が解析結果に影響していない. 図 8 第 1 近似法とフーリエ変換による解析の比較 5. まとめ 本研究では, フーリエ変換による緩和スペクトルの 解析手法を示し, その妥当性を数値実験により検証し た. 以下に結論を述べる. ・ 周波数領域においての誤差を最大限排除すること で解析手法の有効性が確認された. ・ 解析過程でデータに補助区間を設けることで, 正 確な緩和スペクトルが解析できることを確認した. ・ ノイズのある静的試験データに対しても高い精度 で緩和スペクトルが解析できる可能性を示唆した. 参考文献
1) Ferry J D : Viscoelastic Properties of Polymers, 1970 2) Michel Renardy : On the use of Laplace transform
inversion for reconstruction of relaxation spectra, 2008 3) 大崎順彦 : 新・地震動のスペクトル解析入門, 1994 4) 桑村 仁 : 木材の応力緩和試験における技術的問題点と対
策, 日本建築学会構造論文集, 2012
5) Alfrey, T. and Doty, P. : The Methods of Specifying the Properties of Viscoelastic Materials,1945
-1 3 -5 5 15 case 2 0 10 104 1 10-4 10-8 10-12 108 case 2 vh= 4.28 vh= 1.26 -1 20 -5 0 5 10 15 20 case 4 0 10 104 1 10-4 10-8 10-12 108 case 4 vh= 4.96 vh= 0.75
: Theoretical Values : With Auxiliary Section : Without Auxiliary Section
0 5 10 -5 5 15 0 5 10 -5 5 15
: Theoretical Values : Regular procedure : Approximation by Alfrey and Doty
