非黏滯阻尼系統頻響函數靈敏度分析的改進計算方法

史俊磊， 張 磊， 丁 喆

(1. 武漢科技大學 冶金裝備及其控制教育部重點實驗室， 武漢 430081；2. 武漢科技大學 機械傳動與制造工程湖北省重點實驗室， 武漢 430081)

靈敏度分析可以定量地確定系統參數改變對系統動態特性的影響，在結構優化、模型修正和參數識別等領域有著廣泛應用[1-2]。黏性阻尼模型由于其模型簡單和計算方便的優點已被廣泛研究和應用[3]，但隨著各類智能材料和高阻尼復合材料的發展，原有的黏性阻尼模型假設已不能準確地描述該類結構系統的阻尼特性。卷積型非黏滯阻尼模型由于假設阻尼力與質點速度的時間歷程相關，其數學表達式為阻尼核函數與質點速度的卷積,因而能夠準確地描述黏彈性阻尼材料的遲滯效應。目前，非黏滯阻尼系統的研究主要針對其特征值分析、時域響應分析和參數識別等[4-7]方面，其靈敏度分析也主要集中在特征靈敏度分析和時域響應靈敏度分析[8-9]，而非黏滯阻尼系統頻響函數(frequency response function, FRF)靈敏度的研究相對較為缺乏。

FRF靈敏度分析可以評估系統參數對FRF的影響程度。一般來說，結構系統靈敏度的計算方法分為三類：有限差分法[10]、直接微分法和伴隨變量法[11]。其中，有限差分法理論簡單，但在某些情況下計算效率較低。當設計變量數目大于約束數目時，伴隨變量法具有更高的計算效率，反之則直接微分法的計算效率較高。針對FRF靈敏度分析，Lin[12]提出了一種加權FRF靈敏度法，該方法即使在模型誤差較大的情況下，仍具有良好的收斂性能。Qu[13]基于冪級數展開和模態疊加法提出無阻尼系統FRF的模態加速法以及FRF靈敏度的雙模態加速法，減小該系統FRF及其靈敏度的模態截斷誤差。Wu等[14]基于Sturm定理得到參與FRF靈敏度計算的模態數目，并利用激勵頻率收斂冪級數的部分和近似被截斷高階模態影響，進而推導出無阻尼系統的FRF靈敏度表達式。Wu等[15]根據Sherman-Morrison公式，提出了一種僅利用初始FRF和參數攝動來計算靈敏度的方法，該方法避免了重復有限元分析，提高了計算效率。Erik等[16]針對黏滯阻尼系統，提出了一種高計算效率和高精度計算傳遞函數靈敏度的方法。Shadan等[17]利用靈敏度方法，將FRF變化與結構參數變化關聯，提出了一個靈敏度方程來減小FRF數據不完全性的不利影響。Xiao等[18]針對非經典阻尼系統，在僅使用感興趣頻率范圍內模態的情況下推導出FRF靈敏度的表達式，該方法在不涉及狀態空間方程的基礎上修正模態截斷誤差。Zhu等[19]考慮多重擾動對系統動剛度矩陣的影響，采用Sherman-Morrison-Woodbury和有限差分法對系統FRF靈敏度進行了求解。上述FRF靈敏度求解方法均建立在無阻尼或黏性阻尼模型的基礎上，無法直接應用于非黏滯阻尼模型系統。因此亟需提出一種適用于非黏滯阻尼系統FRF靈敏度的計算方法。

本文研究具有卷積型非黏滯阻尼模型結構系統FRF靈敏度問題。非黏滯阻尼系統的特征值求解復雜，尤其對于大型結構系統，特征值分析過程非常耗時。為提高計算效率，通常只利用有限數目的低階模態近似求解其動力響應和靈敏度。但忽略高階模態不可避免地引入模態截斷誤差，對非黏滯阻尼系統動力響應和靈敏度的計算精度影響尤為顯著。因此必須提出一種適用于非黏滯阻尼結構系統FRF靈敏度求解的改進方法，使其能夠補償高階截斷模態的影響，進而提高系統FRF靈敏度計算的精度。解決該問題的難點和關鍵在于建立非黏滯阻尼系統模態和系統矩陣間的關系。為消除非黏滯阻尼系統FRF靈敏度的模態截斷誤差，本文利用Neumann級數展開定理建立結構系統模態與系統矩陣間的關系，考慮其第一項和前兩項近似，將系統高階模態對FRF靈敏度的貢獻用系統矩陣和低階模態進行表達，推導非黏滯阻尼系統FRF靈敏度一階近似法和二階近似方法的表達式。通過數值算例驗證了所提出方法的有效性和準確性。結果表明，與傳統方法相比，本文所提出的改進計算方法兼顧了計算精度和效率，更適用于非黏滯阻尼結構系統FRF靈敏度的求解。

1 理論背景

卷積型非黏滯阻尼模型的阻尼力與質點速度的時間歷程相關，數學表達式為質點速度與某一阻尼核函數的卷積

(1)

g(t)=Cnvg(t)

(2)

式中：t為時間；τ為遲滯時間變量；Fd(t)∈N為阻尼力向量；N為速度向量；Cnv∈N×N為非黏滯阻尼系數矩陣；N為系統的自由度；g(t)為阻尼核函數。式(1)所表示的卷積型非黏滯阻尼模型適用于線性系統，它不僅考慮了當前速度的影響，還考慮了非局部時間的速度影響，更加準確地描述了黏彈性阻尼材料的遲滯效應[20]。

理論上，任何使系統能量耗散函數非負的表達式都可作為卷積型非黏滯阻尼模型的核函數。因此阻尼核函數不唯一，通過選取不同阻尼核函數可適應不同黏彈性材料的阻尼特性，常用的阻尼核函數模型有：Biot阻尼模型[21]、指數阻尼模型和Golla-Hughes-McTavish(GHM)阻尼模型等。特別地，當選取的核函數滿足g(t-τ)=Cδ(t-τ)(C∈N×N是黏性阻尼矩陣，δ(t)為狄拉克函數)時，卷積型非黏滯阻尼模型將退化為傳統的黏性阻尼模型。因此，黏性阻尼模型可看作該非黏滯阻尼模型的一種特殊形式。

含卷積型非黏滯阻尼的N自由度系統時域運動方程可表示為

(3)

其初始條件為

(4)

式中：M和K∈N×N分別是質量矩陣和剛度矩陣；f(t)∈N是外激勵向量；x(t)和N分別是位移和加速度向量。

對式(3)進行拉普拉斯變換得

D(s)X(s)=F(s)

(5)

(6)

式中：λj和φj分別為第j階特征值和特征向量；m為非黏滯阻尼系統的特征值數目m由兩部分組成

m=2N+P，P>0

(7)

式中：2N為彈性模態數目，彈性模態是以復共軛對的形式出現；P為非黏性模態數目，非黏性模態的形式為實數。由式(7)可知，雖然系統有N自由度，但其特征解的個數大于2N。這是非黏滯阻尼系統和黏性阻尼系統的一個重要區別，也是導致非黏滯阻尼系統特征值求解耗時的原因之一。

2 非黏滯阻尼系統FRF靈敏度的求解方法

本章將根據上述非黏滯阻尼系統的運動方程，推導出傳統求解非黏滯阻尼系統FRF靈敏度的表達式。

2.1 直接頻響法

由式(5)可知

X(s)=D-1(s)F(s)=H(s)F(s)

(8)

其中，H(s)為系統FRF，它可表示為

H(s)=D-1(s)=(s2M+sG(s)+K)-1

(9)

該方法稱為直接頻響法(direct frequency response method, DFRM)，將式(9)關于設計變量p求導得到系統FRF靈敏度的表達式

(10)

2.2 模態疊加法

傳統求解系統FRF靈敏度的方法除了DFRM外還有模態疊加法(modal superposition method，MSM)。MSM是在系統特征方程求出系統模態的基礎上，利用留數定理推導出系統FRF矩陣H(s)。當系統全部模態參與計算時，MSM可以獲得精確解，但其計算效率較低，并且在某些情況下，無法求出系統全部模態。因此，通常只考慮低階模態，忽略高階模態對FRF的影響

(11)

(12)

H(s)H(s)-1=I

(13)

式中，I為單位矩陣。同時對式(13)兩邊關于設計變量p求導，并考慮H(s)=D(s)-1得到

(14)

將式(11)代入式(14)得到系統FRF靈敏度表達式

(15)

由式(15)可知，MSM計算出FRF靈敏度的表達式中存在兩個模態截斷項。因此，模態截斷對系統FRF靈敏度分析的影響較大。由模態截斷引起的模態截斷誤差可以定量地表示為

(16)

當L?m時，MSM計算出的FRF靈敏度表達式不僅不準確，且在某些特殊情況下，計算結果是錯誤的。為保證系統FRF靈敏度的計算精度，MSM需要保留較高階模態，但相應的計算效率降低。因此，需要提出一種兼顧計算精度和效率的非黏滯阻尼系統FRF靈敏度的計算方法。

3 非黏滯阻尼系統FRF靈敏度的改進計算方法

在本章將提出兩種考慮高階模態影響的FRF靈敏度的改進計算方法。

3.1 非黏滯阻尼系統模態與系統矩陣關系

本小節介紹一種利用Neumann級數展開定理建立的非黏滯阻尼系統模態與系統矩陣間的關系[22]。

Neumann級數展開定理的表達式為

(I+A)-1=I-A+A2-A3+A4-…

(17)

將式(11)表示為矩陣和向量相乘的形式，并考慮所有階模態得

H(s)=-UΘ-1Λ-1(I-sΛ-1)-1UT

(18)

其中I∈m×m，Λ=diag[λ1,λ2,…,λm]，U=[φ1,φ2,…,φm]和Θ=diag[θ1,θ2,…,θm]。利用Neumann級數展開定理將式(I-sΛ-1)-1展開得

(19)

將式(19)代入式(18)，可以得到

H(s)=-UΘ-1Λ-1UT-sUΘ-1Λ-2UT-

s2UΘ-1Λ-3UT-s3UΘ-1Λ-4UT-…=

(20)

由式(9)可知，FRF還可表示為

H(s)=D-1(s)=(s2M+sG(s)+K)-1=

[I+sK-1(G(s)+sM)]-1K-1

(21)

其中I∈N×N，根據Neumann級數展開定理將[I+sK-1(G(s)+sM)]-1展開并代入到式(21)可以得到

H(s)=K-1-sK-1(G(s)+sM)K-1+[sK-1(G(s)+sM)]2K-1-[sK-1(G(s)+sM)]3K-1+…

(22)

將式(22)展開，并按s的次數由低到高排列可以得到

H(s)=K-1-sK-1G(s)K-1+s2K-1G(s)K-1G(s)K-1-s2K-1MK-1+…=K-1-sK-1G(s)K-1+s2K-1[G(s)K-1G(s)-M]K-1+…=

(23)

其中：

(24)

因此，非黏滯阻尼系統矩陣與模態間關系

-UΘ-1Λ-k-1UT=Γk

(25)

由式(25)可知，由于阻尼核函數G(s)與頻率相關，當k≥2時，所有展開項都會受前項的影響，故其不能再進一步近似等效。因此，考慮非黏滯阻尼系統矩陣與模態間的前兩項近似，并引出其關系

-UΘ-1Λ-1UT=K-1

(26)

-UΘ-1Λ-1UT-UΘ-1Λ-2UT=K-1-K-1G0K-1

(27)

3.2 FRF靈敏度改進計算方法

在上述關系的基礎上，本小節將推導考慮高階模態影響的系統FRF靈敏度的表達式。將式(18)分為可用的低階模態和被截斷的高階模態兩部分

(28)

式中，H表示被截斷的高階模態數目。利用Neumann級數展開將式(28)右側的第一項展開并取前兩階近似得

(29)

聯立式(26)～(29)將被截斷的高階模態用系統低階模態和系統矩陣表示為

(30)

由于高階模態對系統FRF的影響是利用非黏滯阻尼系統模態與系統矩陣關系的前兩項近似表達，因此稱該方法為二階近似法(second-order approximation method, SAM)。聯立式(11)和(30)得到考慮高階模態影響的系統FRF表達式

(31)

通過合并同類項，式(31)可進一步化簡為

(32)

將式(32)代入式(14)得到利用SAM方法計算的非黏滯阻尼系統FRF靈敏度表達式

(33)

SAM計算得到頻響函數對應的模態截斷誤差為

(34)

其中，下標H-SAM表示利用SAM計算出系統頻響函數，其頻響函數靈敏度的模態截斷誤差為

(35)

其中，下標Hp-SAM表示利用SAM計算出系統頻響函數關于p的靈敏度，如果用非黏滯阻尼系統模態與系統矩陣關系第一項近似FRF高階模態的貢獻量可以得到

(36)

聯立式(11)和(34)并化簡可以得到一階近似法(first-order approximation method，FAM)的表達式

(37)

將式(35)代入式(14)得到由FAM計算出的非黏滯阻尼系統FRF靈敏度表達式

(38)

FAM計算得到頻響函數對應的模態截斷誤差為

(39)

其中，下標H-FAM表示利用FAM計算出系統頻響函數，其頻響函數靈敏度的模態截斷誤差為

(40)

其中，下標Hp-FAM表示利用FAM計算出系統頻響函數關于p的靈敏度，與MSM相比，FAM和SAM方法考慮高階模態對非黏滯阻尼系統FRF靈敏度的影響，但計算更加復雜。為了更清晰地表達FAM和SAM的計算流程并便于編程，其算法流程圖如圖1所示。

圖1 FAM和SAM算法流程圖Fig.1 Algorithm flow chart of FAM and SAM

4 例子與結論

4.1 質量-彈簧系統

在本節，通過如圖2所示的含非黏滯阻尼模型的N自由度系統(N可以根據需要進行調整)驗證本文提出的非黏滯阻尼系統FRF靈敏度的FAM和SAM的有效性，并將其與DFRM和MSM進行對比，討論四種方法的計算精度和效率。

本例選取雙指數阻尼核函數，含該阻尼模型的N自由度運動方程表達式為

(41)

式中：c為非黏滯阻尼參數；μ1和μ2均為松弛參數。

圖2 含非黏滯阻尼器的N自由度系統Fig.2 N DOF non-viscously damped system

4.1.1 計算精度

本小節選取自由度數目N=4，系統參數：m1=100 kg，m2=200 kg，m3=300 kg,m4=400 kg和ke=50 000 N/m。非黏滯阻尼模型的參數為：c=50 Ns/m，μ1=50 s-1和μ2=70 s-1。

通過求解狀態空間方程的特征方程可以得到系統的特征值，該系統的特征值如表1所示。可以注意到，本例中的特征值個數為16，其中彈性模態的個數為8，非黏性模態的個數P=r×q=4×2=8(其中r為阻尼系數矩陣的秩，q為頻域內阻尼核函數分母的最高次冪)。由表1可知，彈性模態的形式為復數，其實部表示系統的衰減系數，虛部表示系統的阻尼振動頻率。由于彈性模態的實部與非黏性模態均是負數，所以該系統是穩定的。

表1 四自由度非黏滯阻尼系統的特征解

為了驗證本文提出方法的有效性和準確性，首先選取非黏滯阻尼參數c為設計變量。分別利用DFRM，MSM，FAM和SAM計算該四自由度非黏滯阻尼系統在0～40 rad/s范圍內FRF靈敏度。其中，將DFRM計算得到的FRF靈敏度的結果視為參考值。為了方便起見，圖中只展示系統FRF靈敏度矩陣中第一對角元素H11。H11表示輸出和輸入均為1號自由度時得到的系統FRF。本節中各方法所使用的頻率步長均為0.2 rad/s。

圖3為三種方法選取全部彈性模態計算的FRF靈敏度相對誤差的結果。其中，DFRM計算的實際頻響函數靈敏度也展示在圖3中。本文中，頻響函數靈敏度相對誤差ER的計算公式為

(42)

從圖3可以看出，MSM，FAM和SAM計算出的FRF靈敏度在選取頻率范圍內誤差均在允許范圍內。FAM和SAM的振幅相對誤差明顯小于MSM的振幅相對誤差。SAM的相位相對誤差明顯小于FAM的相位相對誤差。因此，在所比較的三種FRF靈敏度計算方法中，SAM計算精度最高。

圖3 全部彈性模態近似頻響函數靈敏度的相對誤差Fig.3 The relative error of FRF sensitivity matrixapproximated in terms of all elastic modes.

為分析模態截斷誤差對三種方法計算精度的影響，選取第一對(前兩階)彈性模態參與計算FRF的靈敏度，這就意味著系統的模態從第三階開始被截斷，其結果如圖4所示。可以看出，當頻率小于10 rad/s時，MSM的計算結果與參考值仍存在較明顯的誤差，而FAM和SAM的計算結果相較更接近參考值。當頻率大于10 rad/s時，三種方法的FRF靈敏度結果均與參考值存在顯著差異。

圖4 第一對彈性模態近似頻響函數靈敏度Fig.4 The sensitivity of FRF matrix approximated in terms ofthe first pair elastic modes.

圖5是三種方法選取前兩對彈性模態計算出的FRF靈敏度的結果。與圖4結果類似，FAM和SAM的計算精度仍高于相同情形下MSM得到的結果。由于選取的最高階模態頻率為16.71 rad/s，因此FAM和SAM的計算結果在該頻率以內的計算結果與參考值一致，而MSM計算的靈敏度存在較為明顯的誤差。

圖5 前兩對彈性模態近似頻響函數靈敏度Fig.5 The sensitivity of FRF matrix approximated in terms ofthe first two pair elastic modes.

進一步增加參與計算FRF靈敏度的模態數目，考慮前三對彈性模態的各方法計算結果如圖6所示。可以看出，相較于僅考慮第一對和前兩對的靈敏度計算結果，三種方法的FRF靈敏度計算結果均有所提高。通過對比圖4～圖6可以看出，當利用本文提出的FAM和SAM計算非黏滯阻尼系統FRF靈敏度時，若希望得到一定頻率范圍內的準確結果，則需要該頻率范圍內的所有彈性模態參與計算。

圖6 前三對彈性模態近似頻響函數靈敏度Fig.6 The sensitivity of FRF matrix approximated in terms ofthe first three pair elastic modes

為了更加直觀地展示三種方法的準確性，圖7為三種方法選取前三對彈性模態計算出FRF靈敏度的相對誤差。可以看出，各方法的幅值相對誤差高于相應的相位誤差。當選取相同階模態時，FAM和SAM的FRF靈敏度計算結果幾乎相同，但均明顯高于MSM的結果。圖8給出SAM選取不同階模態計算出FRF靈敏度相對誤差。可以看出，增加參與計算的模態數目可大幅降低其FRF靈敏度的相對誤差。結果表明，相較于傳統的MSM方法，由于本文所提出的FAM和SAM方法考慮了非黏滯阻尼系統高階模態截斷誤差的影響，因此能夠大幅提高非黏滯阻尼系統FRF靈敏度的求解精度。

圖7 選取前三對彈性模態計算出頻響函數靈敏度的相對誤差

圖8 SAM使用不同階彈性模態計算FRF靈敏度相對誤差

4.1.2 計算效率

為了比較DFRM，MSM，FAM和SAM的計算效率， 本小節通過圖2所示的N自由度非黏滯阻尼系統，比較四種方法求解該系統頻響函數H11靈敏度所消耗的時間。假定頻率范圍為0～5 rad/s，頻率步長為0.1 rad/s，選取系統自由度N的變化范圍0～1 500。

圖9表示使用不同計算方法的FRF靈敏度計算時間隨自由度增加的變化曲線。可以看出，當自由度小于600時，四種方法的計算時間相差不大。但隨著自由度數目繼續增加，DFRM的計算時間增加迅速，而FAM和SAM的計算時間增長相對平穩。這是因為在每個外激勵頻率下DFRM要對系統動剛度矩陣求逆，而系統自由度增加必然會引起系統動剛度矩陣維數的增加，系統動剛度矩陣求逆消耗的時間就越長。而FAM和SAM避免了矩陣求逆，僅在MSM的基礎上考慮高階模態對系統頻響函數靈敏度的影響。因此，對于大自由度系統，FAM和SAM相較于DFRM在計算效率上具有明顯的優勢，且FAM的效率稍高于SAM。這是因為FAM和SAM分別考慮利用系統矩陣與模態關系的第一項和前兩項近似高階模態對系統頻響函數靈敏度的影響。

圖9 四種方法計算時間隨自由度數目變化圖Fig.9 Computational time of different DOF with four methods

4.2 軸向振動桿

本節考慮一個具有雙指數阻尼的固定-自由軸向振動桿。如圖10所示，該軸向振動桿可根據需求劃分為N+1個單元。軸向振動桿單元的單元剛度矩陣和單元質量矩陣如下

(43)

式中：ρ表示密度；A為桿的橫截面積；E表示材料的楊氏模量；le=L/N為單元桿的長度。各參數分別為A=6.25×10-4m2，E=2.1×1011N/m2,ρ=7.8×103kg/m3，桿的總長L=4 m。整體質量矩陣M和整體剛度矩陣K可以根據經典的有限元方法組裝得到。

假設軸向振動桿含雙指數阻尼模型的表達式為

Cnv=μ1e-μ1(t)Cnv1+μ2e-μ2(t)Cnv2

(44)

其中，松弛因子μ1和μ2為

(45)

Cnv1=αM,Cnv2=βK

(46)

其中

(47)

圖10 含雙指數阻尼模型的軸向振動桿示意圖Fig.10 Axially vibrating rod with double-exponentialdamping model

4.2.1 計算精度

為了討論不同方法的計算精度，本節將該軸向振動桿劃分為80個單元，即N=80。分別使用DFRM、MSM、FAM和SAM計算該軸向振動桿自由端的頻響函數H11關于ρ的靈敏度分析，選取的頻率范圍為：0～10 000 rad/s，選取的頻率步長為10 rad/s。本算例仍然以DFRM計算出的頻響函數靈敏度為參考值。

圖11給出使用前兩對彈性模態計算出系統頻響函數靈敏度，各方法相對于DFRM的相對誤差如圖12所示。在所選頻率范圍內，FAM和SAM計算出頻響函數靈敏度的精度優于相同情況下MSM的計算精度。對于FAM和SAM，二者幅值的相對誤差相差無幾，而SAM的相位相對誤差明顯比FAM的相位相對誤差小。因此，就計算精度而言，SAM的計算精度高于FAM計算精度。這與質量-彈簧系統算例得到的結論相同。

圖11 前兩對彈性模態近似頻響函數靈敏度

4.2.2 計算效率

本小節通過改變軸向振動桿的自由度討論DFRM、MSM、FAM和SAM的計算效率。將軸向振動桿自由度劃分為N=0～800，自由度間隔為80。假定頻率范圍為0～10 000 rad/s，頻率步長為100 rad/s。圖13是四種方法計算效率對比圖。可以看出，FAM和SAM的計算效率明顯優于DFRM，而FAM的計算效率高于SAM。各方法的計算效率排序與質量-彈簧系統算例得到的結論相同。

圖12 前兩對彈性模態近似頻響函數靈敏度的相對誤差Fig.12 The relative error of FRF sensitivity approximated in terms

圖13 不同方法計算時間隨自由度數目變化圖Fig.13 Computational time of different DOF with different methods

5 結 論

本文提出了具有卷積型非黏滯阻尼模型結構系統頻響函數靈敏度的改進計算方法。利用復模態疊加法和直接微分法推導出了非黏滯阻尼系統頻響函數靈敏度表達式。但若只用有限數目的低階模態求解其頻響函數靈敏度會引入模態截斷誤差，在某些情形下會嚴重影響非黏滯阻尼系統頻響函數靈敏度的求解精度。解決該問題的關鍵在于建立的系統模態和系統矩陣間的關系。本文利用Neumann級數展開定理，將系統高階模態對系統頻響函數靈敏度的貢獻用系統矩陣和低階模態進行表達，進而推導出了兩種考慮高階模態影響的系統頻響函數靈敏度的改進計算方法：一階近似法和二階近似法。通過數值算例驗證兩種方法的有效性和準確性，結果表明：

(1) 一階近似法和二階近似法能得到滿意的系統頻響函數靈敏度計算結果，因此本文提出的兩種改進計算方法可有效地求解非黏滯阻尼系統頻響函數靈敏度。

(2) 一階近似法和二階近似法由于考慮了高階模態對系統頻響函數靈敏度的影響，因此其計算精度優于模態疊加法；從計算效率角度看，改進的一階近似法和二階近似法避免了矩陣求逆，因此其計算效率高于直接頻響法。綜上，本文提出的兩種改進計算方法兼顧了計算精度和效率，是求解多自由度非黏滯阻尼系統頻響函數靈敏度問題的理想選擇。

(3) 由于一階近似法和二階近似法分別考慮利用系統矩陣與模態關系的第一項和前兩項近似高階模態對系統頻響函數靈敏度的影響，因此二階近似法比一階近似法的計算精度更高，而計算效率卻低于一階近似法。綜合考慮計算精度和效率，建議根據需要選取不同計算方法：當需要較高計算精度時，選取二階近似法；當需要較高計算效率時，選取一階近似法。