知識科普 - 結構物風致振動之頻率域解析理論

參考自「風工程理論與應用第四章結構風致振動之頻率域解析劉明怡教授主筆」

風速、風力與結構物振動反應之關係

風引致結構物的振動行為可分為順風向反應（along-wind response）、橫風向反應（across-wind response）與扭轉向反應（torsional response）等三個分量，其振動機制均有所差異，但均須考量其同時作用於結構物的綜合反應。順風向反應的主因為抖振（buffeting），亦即來流對結構物產生的擾動風力所激發的結構物振動，其方向平行於順風向（Balendra 1993; Dyrbye and Hansen 1997; Holmes 2015; Liu 1991; Simiu 2011; Simiu and Scanlan 1996; Tamura and Kareem 2013; 朱 2006）。橫風向反應與扭轉向反應的主因則是來流在結構物周圍所發生分離、再接觸、尾跡等鈍體空氣動力學現象的結果，多以風洞試驗所迴歸的經驗公式來評估其反應大小。以下針對順風向反應的評估方法進行說明。

紊流(turbulence)風場中之風速、風力與結構振動反應等三個參數隨時間的變化均具有隨機性。因此，定值法（deterministic approach）無法合理地模擬這些現象。為了改善定值法之侷限性，Davenport (1963)從機率法（probabilistic approach）的觀點出發，提出以隨機振動（random vibration）和頻譜分析（spectral analysis）的理論為基礎，執行頻率域分析（frequency domain analysis）的流程。此流程具有以下幾點假設：

(1) 風速、風力和結構物的振動反應的歷時序列（time history）均為定常隨機過程（stationary random process），亦即此三個物理量參數之平均值（mean value）、標準差（standard deviation value）與尖峰值（peak value）等統計參數，皆不隨時間區間的不同而變化（Bendat and Piersol 2010; Newland 1993）。
(2) 假設結構物屬於線彈性系統（linear elastic system），其質量（mass）、阻尼（damping）和勁度（stiffness）等動力參數不隨著時間而改變。
(3) 當結構物的振動行為發生時，材料處於微小變形狀態（Chopra 2012）。

如圖4-1的分析流程所示，係針對順風向反應的標準差進行評估。根據隨機振動學理論，各參數之頻譜密度函數（spectral density function）所圍區域面積即為該參數的變異數（variance），並分別以氣動力導納函數（aerodynamic admittance function）與機械導納函數（mechanical admittance function）連接這些參數之頻譜密度函數。

圖4-1 順風向風力等值靜載重之分析流程

單自由度系統順風向反應評估流程

單自由度系統（single-degree-of-freedom system）如圖4-2所示，其中$m$、$c$ 與 $k$ 分別為系統之質量、阻尼和勁度，$F(t)$ 為順風向風力歷時，$x(t)$ 為順風向位移歷時，$t$ 則為時間，其運動方程式（equation of motion）為：

$$m\ddot{x}(t)+c\dot{x}(t)+kx(t)=F(t)\tag{4-1}$$

式(4-1)可改寫為：

$$\ddot{x}(t)+2\xi_0(2{\pi}n_0)\dot{x}(t)+(2{\pi}n_0)^2x(t)=\frac{F(t)}{m}\tag{4-2}$$

其中

$$n_0=\frac{1}{2\pi}\sqrt{\frac{k}{m}}\tag{4-3}$$

$$\xi_0=\frac{c}{2\sqrt{km}}\tag{4-4}$$

$n_0$ 與 $\xi_0$ 則分別為結構物的自然頻率（natural frequency）和阻尼比（damping ratio）。

圖4-2 單自由度系統示意圖

風速歷時 $u(t)$ 為平均風速 $\bar{u}$ 與擾動風速歷時 $u'(t)$ 之疊加，亦即：

$$u(t)=\bar{u}+u'(t)\tag{4-5}$$

在準穩定假設（quasi-steady assumption）的條件下，$F(t)$ 可直接由 $u(t)$ 推導而得，亦即：

$$F(t)=\frac{1}{2}{\rho}BD{C_D}{u^2}(t)\tag{4-6}$$

其中 $\rho$ 與 $C_D$ 分別為空氣密度與順風向風力係數（或阻力係數drag coefficient）。$B$ 與 $D$ 可視為結構物迎風面寬度及高度，故式中 $BD$ 表示為順風向風力作用面積。將式(4-5)代入式(4-6)後，可得：

$$F(t)=\frac{1}{2}{\rho}BD{C_D}\left[\bar{u}^2+2\bar{u}{u^\prime}(t)+{u^\prime}^2(t)\right]\tag{4-7}$$

就風工程應用觀點而言， $u^\prime(t)/\bar{u}$ 甚少超過0.2，故式(4-7)之 ${u^\prime}^2(t)$ 甚小，可忽略不計（Simiu and Scanlan 1996），亦即：

$$F(t)\approx\frac{1}{2}{\rho}BD{C_D}\bar{u}^2+{\rho}BD{C_D}\bar{u}{u^\prime}(t)=\bar{F}+F^\prime(t)\tag{4-8}$$

其中 $\bar{F}$ 和 $F^\prime(t)$ 分別為順風向風力平均值與擾動歷時。式(4-8)顯示 $F(t)$ 為 $\bar{F}$ 和 $F^\prime(t)$ 的疊加。假設 $u^\prime(t)$ 與 $F^\prime(t)$ 在時間區間 $T$ 內，皆為平均值等於0的定常隨機過程，則擾動風力歷時的自相關函數（autocorrelation function）$R_F(\tau)$ 可推導如下：

$$R_F(\tau)=\lim_{T \to \infty}\frac{1}{T}\int_{-T/2}^{T/2} F^\prime(t)F^\prime(t+\tau) dt = \overline{F^\prime(t) F^\prime(t+\tau)} = (\rho BD C_D \bar{u})^2 \overline{u^\prime(t) u^\prime(t+\tau)}\tag{4-9}$$

其中，$\tau$ 為時間間隔，(─)則表示整體平均（ensemble average）。將式(4-9)進行傅利葉轉換(Fourier transform)後，可得

$$S_F(n)=2\int_{-\infty}^{\infty}R_F(\tau) \cos{2\pi n} \tau d\tau = (\rho BD C_D \bar{u})^2 S_u(n)=\frac{4\overline{F}^2}{\overline{u}^2}S_u(n)\tag{4-10}$$

其中， $S_u(n)$ 和 $S_F(n)$ 分別為擾動風速歷時與擾動風力歷時的頻譜密度函數，$n$ 則為頻率。就建築結構物設計觀點而言，可使用高程 $z$ 處之擾動風速歷時的頻譜密度函數 $S_u(z,n)$ 取代式(4-10)中之 $S_u(n)$（Kaimal et al. 1972），亦即：

$$\frac{n S_u (z,n)}{u_*^2}=\frac{200f}{(1+50f)^{5/3}}\tag{4-11}$$

其中，

$$u_*^2=\frac{\sigma_u^2(z)}{\beta}\tag{4-12}$$

$$f=\frac{nz}{\overline{u}(z)}\tag{4-13}$$

$u_\ast$ 和 $f$ 分別為摩擦速度(friction velocity)與無因次頻率，$\bar{u}(z)$ 和 ${\sigma_u}^2(z)$ 分別為高程 $z$ 處的平均風速與擾動風速歷時之均方值，$\beta$ 則為和粗糙長度（roughness length）有關的係數。圖4-3為 $S_u(z,n)$ 與 $n$ 之關係曲線，顯示擾動風速能量集中在低頻區域，且隨著 $n$ 增加而遞減，呈現寬頻（broad bandwidth）特性；此外，$S_u(z,n)$ 所圍區域面積等於 ${\sigma_u}^2(z)$，亦即：

$$\sigma_u^2=\int_{0}^{\infty} S_u(z,n) dn\tag{4-14}$$

圖4-3 擾動風速歷時之頻譜密度函數

紊流風場充滿各種尺度的渦漩（eddy），結構物之存在將干擾這些渦漩，並改變其迎風面上擾動風速的相關性，進而影響渦漩作用於結構物所產生之阻力。因此，必須以氣動力導納函數 $\chi^2(n)$ 修正式(4-10)（Vickery 1965），亦即：

$$S_F(n)=\frac{4\overline{F}^2}{\overline{u}^2}\chi^2(n)S_u(n)\tag{4-15}$$

其中

$$\chi(n)=\frac{1}{1+\left(\frac{2\sqrt{BD}n}{\overline{u}}\right)^{4/3}}\tag{4-16}$$

圖4-4為 $\chi^2(n)$ 和 $n$ 的關係曲線，顯示低頻區域之渦漩具有較長波長（wavelength），可完整地包覆結構物並產生阻力（順風向風力），此時，結構迎風面上擾動風速可視為完全相關，故 $\chi^2(n)$ 趨近於 1；此相關性隨著渦漩頻率增加（渦漩波長減短）而降低，渦漩作用於結構物所產生的阻力亦隨之減小，故 $\chi^2(n)$ 呈遞減趨勢；高頻區域的 $\chi^2(n)$ 則趨近於 0，此時，結構迎風面上擾動風速之相關性幾乎完全消失，導致波長較短的渦漩無法有效地在結構物上產生阻力。

圖4-4 氣動力導納函數

圖4-5則為 $S_F(n)$ 與 $n$ 之關係曲線，顯示擾動阻力能量集中在低頻區域，且隨著 $n$ 增加而遞減，呈現寬頻特性；此外，$S_F(n)$ 所圍區域面積等於擾動風力歷時的均方值 ${\sigma_F}^2$，亦即：

$$\sigma_F^2=\int_{0}^{\infty} S_F(n) dn\tag{4-17}$$

$x(t)$ 為結構物在順風向之平均位移 $\bar{x}$ 和擾動位移歷時 $x^\prime(t)$ 的疊加，亦即：

$$x(t)=\bar{x}+x^\prime(t)\tag{4-18}$$

其中

$$\bar{x}=\frac{\bar{F}}{k}\tag{4-19}$$

圖4-5 擾動風力歷時之頻譜密度函數

假設圖4-2之 $F(t)$ 具有諧和力（harmonic force）的形式，亦即：

$$F(t)=\bar{F}\sin{2\pi nt}\tag{4-20}$$

則式(4-1)之結構物順風向穩態位移函數（steady-state displacement function）$x_{st}(n)$ 可推導如下

$$x_{st}(n)=\frac{\bar{F}}{k}\frac{1}{\sqrt{\left[1-\left(\frac{n}{n_0}\right)^2\right]^2+\left[4\xi_0^2\left(\frac{n}{n_0}\right)\right]^2}}\tag{4-21}$$

式(4-21)與式(4-19)的比值之平方定義為機械導納函數 $\left|H(n)\right|^2$（Chopra 2012），亦即：

$$\left|H(n)\right|^2=\left[\frac{x_{st}(n)}{\overline{x}}\right]^2=\frac{1}{\left[1-\left(\frac{n}{n_0}\right)^2\right]^2+4\xi_0^2\left(\frac{n}{n_0}\right)^2}\tag{4-22}$$

圖4-6為 $\zeta_0$ 介於 0 和 1 之間時，$\left|H(n)\right|^2$ 與 $n$ 的關係曲線，顯示低頻區域之各個 $\zeta_0$ 對應的 $\left|H(n)\right|^2$ 趨近於 1，亦即 $x_{st}(n)$ 接近 $\bar{x}$；高頻區域之各個 $\zeta_0$ 對應的 $\left|H(n)\right|^2$ 則趨近於 0，亦即 $x_{st}(n)$ 接近 0；當結構物之振動頻率接近自然頻率（$n$ 趨近於 $n_0$）時，系統發生共振（resonance），在大多數土木材料的 $\zeta_0$ 均小於 0.05 之前提下，具有明顯尖峰的窄頻（narrow bandwidth）特性，亦即 $x_{st}(n)$ 遠大於 $\bar{x}$，此現象隨著 $\zeta_0$ 遞減而更趨顯著。

圖4-6 機械導納函數

將式(4-15)和式(4-22)合併後，可推導結構順風向擾動位移歷時之頻譜密度函數 $S_x(n)$，亦即：

$$S_x(n)=\frac{1}{k^2}\left|H(n)\right|^2S_F(n)\tag{4-23}$$

圖4-7為 $S_x(n)$ 與 $n$ 的關係曲線，顯示在 $\zeta_0$ 小於 0.05 之前提下，$S_x(n)$ 所圍區域包含背景反應（background response）區和共振反應（resonant response）區兩部分。前者的結構順風向擾動位移能量集中在低頻區域，且隨著 $n$ 增加而遞減，呈現寬頻特性，與 $S_F(n)$ 之趨勢相同，因此，該區能量來自於紊流產生的擾動阻力；後者之結構順風向擾動位移能量則集中在 $n$ 趨近於 $n_0$ 處，具有明顯尖峰的窄頻特性，和 $\left|H(n)\right|^2$ 之趨勢相同，因此，該區能量來自於結構共振。此外，圖4-8中的 $S_x(n)$ 所圍區域面積等於結構順風向擾動位移歷時之均方值 ${\sigma_x}^2$，亦即：

$$\sigma_x^2=\int_{0}^{\infty} S_x(n) dn = \int_{0}^{\infty} \left|H(n)\right|^2 S_F(n) dn\tag{4-24}$$

圖4-7 結構順風向擾動位移歷時之頻譜密度函數

為了簡化式(4-24)複雜的積分運算，假設背景反應區之各個 $n$ 對應的 $\left|H(n)\right|^2$ 等於 1，且共振反應區之 $S_F(n_0)$ 為常數，則式(4-24)近似於背景反應區與共振反應區兩部分的面積總和，亦即：

$$\sigma_x^2 \approx \frac{1}{k^2} \int_{0}^{\infty} S_F(n) dn + \frac{S_F(n_0)}{k^2} \int_{0}^{\infty} \left|H(n)\right|^2 dn = \frac{1}{k^2} \int_{0}^{\infty} S_F(n) dn + \frac{S_F(n_0)}{k^2} \frac{\pi n_0}{4\xi_0} = \sigma_B^2 + \sigma_R^2 \tag{4-25}$$

其中

$$\int_{0}^{\infty}\left|H(n)\right|^2 dn=\frac{\pi n_0}{4\xi_0}\tag{4-26}$$

可由隨機振動理論推導而得（Crandall and Mark 1963; Newland 1993），${\sigma_B}^2$ 和 ${\sigma_R}^2$ 則分別為結構順風向擾動位移歷時之背景反應均方值與共振反應均方值，兩者分別等於背景反應區和共振反應區的面積。以圖4-1之風引致結構順風向反應的分析流程為基礎，Davenport (1964)提出結構順風向尖峰位移 $\hat{x}$ 之計算式，亦即：

$$\hat{x}=\overline{x}+\hat{x}^\prime+g\sigma_x\tag{4-27}$$

其中，$\hat{x}\prime$ 與 $g$ 分別為結構順風向擾動位移歷時的尖峰值和尖峰因子（peak factor）。假設 $x\prime(t)$ 在 $T$（通常取3,600秒）之間為平均值等於 0 的高斯隨機過程（Gaussian random process），則 $g$ 可推導如下：

$$g=\sqrt{2\ln{\nu T}}+\frac{0.577}{\sqrt{2\ln{\nu T}}}\tag{4-28}$$

其中，$\nu$ 為 $x\prime(t)$ 通過時間座標軸之斜率為正值的跨零點頻率（zero-crossing frequency），在 $\zeta_0$ 小於 0.05 之前提下，$\nu$ 近似於 $n_0$。