黃江濤，高正紅，白俊強(qiáng)，周 鑄，趙 軻

（1.西北工業(yè)大學(xué) 翼型葉柵空氣動(dòng)力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室，西安 710072；2.中國空氣動(dòng)力研究與發(fā)展中心，綿陽 621000）

0 引 言

飛機(jī)氣動(dòng)彈性特性研究傳統(tǒng)上劃為靜氣動(dòng)彈性和顫振特性研究兩個(gè)獨(dú)立部分［1－2］。靜氣動(dòng)彈性特性計(jì)算求解過程往往采用靜變形和定常氣動(dòng)力的交替迭代計(jì)算，最后得出發(fā)散速度或特定飛行條件下的機(jī)翼靜變形。傳統(tǒng)的飛行器設(shè)計(jì)主要是針對(duì)型架外形設(shè)計(jì)。然而在飛行器經(jīng)過氣動(dòng)彈性變形之后，巡航狀態(tài)氣動(dòng)性能與設(shè)計(jì)性能存在一定的區(qū)別。而對(duì)于民用客機(jī)來講首先進(jìn)行巡航外形設(shè)計(jì)，通過靜氣動(dòng)彈性計(jì)算、修正得到型架外形，從而保證型架外形在巡航狀態(tài)飛行時(shí)，在氣動(dòng)彈性作用下氣動(dòng)外形與性能能夠恢復(fù)到巡航設(shè)計(jì)狀態(tài)。在型架修正及分析計(jì)算方面國內(nèi)外采用不同的方法進(jìn)行了不同程度的研究。例如國內(nèi)楊國偉采用稀疏背景網(wǎng)格進(jìn)行動(dòng)網(wǎng)格變形、RBF進(jìn)行流固耦合，從而實(shí)現(xiàn)型架修正以及靜氣彈分析研究［3］，詹浩等采用反設(shè)計(jì)以及簡化結(jié)構(gòu)模型方法進(jìn)行了機(jī)翼型架設(shè)計(jì)及分析計(jì)算［4］，國外C B Allen等采用RBF結(jié)合TFI技術(shù)進(jìn)行了氣動(dòng)彈性的分析計(jì)算［5］，取得了合理的計(jì)算結(jié)果。然而對(duì)于氣動(dòng)彈性計(jì)算來講，變形網(wǎng)格的工作效率一定程度上決定了氣動(dòng)彈性分析計(jì)算效率，尤其在非定常動(dòng)氣動(dòng)彈性計(jì)算中需要不斷的調(diào)用變形網(wǎng)格技術(shù)，這時(shí)網(wǎng)格變形工作效率更為重要。

對(duì)于靜氣動(dòng)彈性來說，其關(guān)鍵技術(shù)為空間網(wǎng)格變形技術(shù)以及CFD／CSD數(shù)據(jù)交換技術(shù)。其中變形網(wǎng)格技術(shù)效率決定著靜氣動(dòng)彈性計(jì)算效率、精度，其穩(wěn)健性決定了能否進(jìn)行飛行器大運(yùn)動(dòng)量下網(wǎng)格變形且能否保持較高網(wǎng)格質(zhì)量，以及所求氣動(dòng)彈性變形外形的光滑程度；CFD／CSD數(shù)據(jù)交換技術(shù)決定著空氣動(dòng)力學(xué)模型與結(jié)構(gòu)模型的數(shù)據(jù)能否進(jìn)行高效合理的交換。

本文基于徑向基函數(shù)與圖映射技術(shù)建立了高效的網(wǎng)格變形技術(shù)，與相關(guān)文獻(xiàn)不同的是，本文在合理提取飛行器表面稀疏點(diǎn)序列前提下，有效的利用多塊網(wǎng)格的塊頂點(diǎn)進(jìn)行Delaunay圖映射單元的構(gòu)建（網(wǎng)格分塊盡量均勻），避免了傳統(tǒng)Delaunay圖單元構(gòu)建時(shí)，相關(guān)點(diǎn)輸入的復(fù)雜操作，并通過RBF技術(shù)對(duì)網(wǎng)格頂點(diǎn)進(jìn)行高效操作，實(shí)現(xiàn)映射單元變形；通過建立共用的“映射曲面”，實(shí)現(xiàn)了CFD／CSD數(shù)據(jù)高效高精度交換，構(gòu)建了穩(wěn)健性強(qiáng)的靜氣動(dòng)彈性計(jì)算平臺(tái)；利用該技術(shù)對(duì)某型支線客機(jī)進(jìn)行了型架外形設(shè)計(jì)、修正，進(jìn)一步通過靜氣動(dòng)彈性分析計(jì)算，驗(yàn)證本文所建立方法的正確性。

1 Delaunay圖映射技術(shù)

本文建立了基于Delaunay圖映射的網(wǎng)格變形技術(shù)。Dirichlet圖是由一組連接兩相鄰點(diǎn)直線的垂直平分線組成的連續(xù)多邊形，而Delaunay三角網(wǎng)是由與相鄰Dirichlet多邊形共享一條邊的相關(guān)點(diǎn)連接而成的三角形，換句話說Dirichlet三角網(wǎng)是Delaunay圖的偶圖。

基于Delaunay圖映射的網(wǎng)格變形技術(shù)基本依據(jù)為：在平面或空間給定包括流場邊界點(diǎn)在內(nèi)的一組點(diǎn)，即可以進(jìn)行唯一的Delaunay圖三角化。在Delaunay圖四面體單元構(gòu)建時(shí)，依據(jù)球面準(zhǔn)則進(jìn)行交換測試形式的優(yōu)化，根據(jù)插入點(diǎn)坐標(biāo)定位點(diǎn)所在格網(wǎng)，如果該格網(wǎng)記錄有三角形編號(hào)，則從該編號(hào)開始在相鄰三角形中定位插入點(diǎn)所在三角形；由近及遠(yuǎn)遍歷格網(wǎng)的相鄰網(wǎng)格，若存在三角編號(hào)則遍歷結(jié)束；在最新生成三角形中定位插入點(diǎn)。從而可以完成計(jì)算區(qū)域的三角化圖覆蓋，進(jìn)一步對(duì)空間計(jì)算域任意點(diǎn)進(jìn)行四面體映射單元定位。

計(jì)算域網(wǎng)格節(jié)點(diǎn)定位完畢后，需建立映射關(guān)系［6］，對(duì)于平面Delaunay三角化背景網(wǎng)格，建立如圖1的計(jì)算網(wǎng)格與平面Delaunay三角單元關(guān)系，對(duì)于任意網(wǎng)格節(jié)點(diǎn)O處于三角形MNQ中，由O和三角形頂點(diǎn)構(gòu)造的三角形MON、QOM以及NOQ的面積分別為S1、S2、S3，令：

則可以建立Delaunay三角單元節(jié)點(diǎn)與計(jì)算網(wǎng)格節(jié)點(diǎn)如下關(guān)系式：

飛行器氣動(dòng)外形變形時(shí)，相當(dāng)于Delaunay圖背景網(wǎng)格三角單元變化，即 ［XMXNXQ］更新為［X′MX′NX′Q］，對(duì)于三維網(wǎng)格變形方法與二維原理一致，如圖2所示。依據(jù)建立的Delaunay圖映射關(guān)系式，計(jì)算網(wǎng)格節(jié)點(diǎn)更新為：

以氣動(dòng)彈性變形后飛行器表面區(qū)域頂點(diǎn)位移量為輸入量，建立RBF徑向基函數(shù)插值系統(tǒng)，進(jìn)行空間區(qū)域頂點(diǎn)位移量插值，操作空間網(wǎng)格區(qū)域頂點(diǎn)運(yùn)動(dòng)，從而進(jìn)行Delaunay圖四面體映射單元的變形，利用式（2）進(jìn)行計(jì)算網(wǎng)格的變形，完成網(wǎng)格更新。而Delaunay映射變換之前首先要完成空間區(qū)域頂點(diǎn)的更新，即采用RBF徑向基函數(shù)技術(shù)，其工作示意圖見圖3。

圖1 網(wǎng)格節(jié)點(diǎn)與Delaunay三角形Fig.1 Grid points and Delaunay triangle

圖2 網(wǎng)格節(jié)點(diǎn)與Delaunay四面體Fig.2 Grid points and Delaunay tetrahedron

圖3 RBF工作示意圖Fig.3 RBF work map

2 RBF徑向基函數(shù)

RBF徑向基函數(shù)基本概念由Buhmann和Wendland提出［7－8］，該技術(shù)在精確插值、神經(jīng)網(wǎng)絡(luò)構(gòu)建、數(shù)據(jù)預(yù)測等方面獲取了廣泛的應(yīng)用［9－13］?；赗BF方法的插值數(shù)學(xué)模型可以表示為：

本文對(duì)p（x）的選取采用多項(xiàng)式方法：

其中M 為徑向基中心矢量的維數(shù)。p（x）的建立不僅能夠精確的描述物體或網(wǎng)格的平移以及旋轉(zhuǎn)運(yùn)動(dòng)方式，且進(jìn)行徑向基插值時(shí)能夠保證輸入點(diǎn)系統(tǒng)“力平衡”及“力矩平衡”［14］。其中，xi為徑向基中心，相當(dāng)于圖3中ci；αi為插值系數(shù)，相當(dāng)于圖3中λi?！瑇－xi‖為x到中心xi的徑向歐氏距離。且徑向基函數(shù)φ具有以下性質(zhì)：

以輸入數(shù)據(jù)為RBF技術(shù)構(gòu)建基礎(chǔ)，可以寫為如下矩陣形式：

其中：Φ ＝［φij］，φij＝φ（‖xi－xj‖），（i，j＝1，2，…，N ）；Φ 為插值矩陣，ΓT＝Φ－1Y，Γ ＝ ［λ1，λ2，…，λN］。從而實(shí)現(xiàn)RBF徑向基函數(shù)框架構(gòu)建，RBF（x）則為通過所有輸入點(diǎn)的連續(xù)可微函數(shù)。

3 CFD／CSD數(shù)據(jù)耦合技術(shù)與靜力學(xué)方程求解

在CFD／CSD數(shù)據(jù)交換時(shí)，本文采用構(gòu)建空氣動(dòng)力模型與結(jié)構(gòu)模型的共用“映射曲面”［15］，根據(jù)機(jī)翼形狀建立一曲面，該曲面由兩個(gè)方向邏輯參數(shù)進(jìn)行描述，網(wǎng)格中第i行，第j列處的節(jié)點(diǎn)坐標(biāo)為（x，y），且該點(diǎn)處的參數(shù)化坐標(biāo)定義為（u ，v），u ＝i，v ＝j(luò)。將映射網(wǎng)格控制點(diǎn)參數(shù)化后得到參數(shù)空間。

對(duì)于CFD／CSD網(wǎng)格交界面上節(jié)點(diǎn)的參數(shù)化就是要計(jì)算該點(diǎn)在參數(shù)空間中的坐標(biāo)，首先利用面積坐標(biāo)方法［15］定位源／目標(biāo)學(xué)科節(jié)點(diǎn)所處映射曲面單元（i，j），進(jìn)一步通過雙線性曲面插值計(jì)算該點(diǎn)在其所在的結(jié)構(gòu)化映射網(wǎng)格單元內(nèi)的位置（i0，j0），則該節(jié)點(diǎn)在參數(shù)空間的位置為（i＋i0，j＋j0）。

本文采用雙三次Bezier曲面［15－16］進(jìn)行源學(xué)科節(jié)點(diǎn)與目標(biāo)學(xué)科節(jié)點(diǎn)變量插值：

其中，U、V為目標(biāo)學(xué)科參數(shù)空間邏輯坐標(biāo)，P為源學(xué)科節(jié)點(diǎn)矩陣，M 為Bezier基函數(shù)系數(shù)矩陣。

圖4 “映射曲面”工作示意圖Fig.4 Mapping surface work map

采用柔度系數(shù)矩陣方法［17－18］求解結(jié)構(gòu)靜力學(xué)方程時(shí)，僅考慮飛機(jī)機(jī)翼的彈性變形，發(fā)動(dòng)機(jī)及其掛架運(yùn)動(dòng)與機(jī)翼連接處保持一致。

其中，i＝（1，2，…，n），常數(shù)cij為柔度影響系數(shù)矩陣。則結(jié)構(gòu)靜力學(xué)方程為：

qs為結(jié)構(gòu)節(jié)點(diǎn)在三個(gè)方向的變形位移矩陣；C為結(jié)構(gòu)柔度矩陣；Fs為作用在節(jié)點(diǎn)上的氣動(dòng)力矩陣。對(duì)于型架構(gòu)型逆計(jì)算來講，只需將式（7）中力方向反向。綜上所述，本文建立的靜氣動(dòng)彈性數(shù)值模擬流程見圖5。

圖5 靜氣動(dòng)彈性計(jì)算流程Fig.5 Static aeroelastic work flow chart

4 某型支線客機(jī)型架外形修正與靜氣彈計(jì)算

以某型支線客機(jī)為型架外形為計(jì)算算例，進(jìn)行型架設(shè)計(jì)后修型，進(jìn)一步進(jìn)行靜氣動(dòng)彈性計(jì)算，驗(yàn)證所求型架構(gòu)型正確性以及本文相關(guān)技術(shù)的效率、可靠性。

該型客機(jī)設(shè)計(jì)巡航狀態(tài)為：

該型客機(jī)結(jié)構(gòu)有限元采用桿梁單元，建立雙梁單塊式結(jié)構(gòu)模型，主要結(jié)構(gòu)元素包括前梁、后梁、翼肋等。利用有限元分析計(jì)算了該結(jié)構(gòu)模型的剛度矩陣?？諝鈩?dòng)力學(xué)模型采用多區(qū)域點(diǎn)對(duì)接結(jié)構(gòu)網(wǎng)格進(jìn)行空間離散，空間流場流場共分為330個(gè)區(qū)域，表面劃分為160個(gè)區(qū)域，共1300萬個(gè)網(wǎng)格單元（半模）。圖6給出了該支線客機(jī)的表面網(wǎng)格分布示意圖。氣動(dòng)分析主控方程為Navier－Stokes方程；空間離散方法為二階Roe格式；湍流模型采用一方程SA模型；采用LU－SGS隱式時(shí)間推進(jìn)。通過線性插值的方法求出設(shè)計(jì)升力系數(shù)下所對(duì)應(yīng)的巡航攻角α＝1.95°。圖7為巡航狀態(tài)下飛機(jī)表面壓力云圖。

根據(jù)CFD數(shù)值模擬結(jié)果，采用文中建立的CFD／CSD共用的“映射曲面”技術(shù)進(jìn)行氣動(dòng)載荷向結(jié)構(gòu)載荷的映射，通過求解靜力學(xué)方程，求出結(jié)構(gòu)節(jié)點(diǎn)位移，再一次通過“映射曲面”技術(shù)將結(jié)構(gòu)位移映射到氣動(dòng)網(wǎng)格上，從而完成型架構(gòu)型的逆向計(jì)算，如圖所示，綠色外形為型架構(gòu)型，淺色為巡航設(shè)計(jì)外形，可以看出機(jī)翼展向下彎較為明顯，表1給出了展向10個(gè)站位型架外形機(jī)翼與巡航外形的相對(duì)扭角。

對(duì)于生產(chǎn)與裝配工藝來講，要求機(jī)翼前緣為一直線且與巡航設(shè)計(jì)外形重合，因此在保持型架構(gòu)型不變的情況下對(duì)機(jī)翼做前緣一致性修正，以滿足裝配工藝要求。

圖6 CFD表面網(wǎng)格Fig.6 CFD surface grid

圖7 巡航狀態(tài)表面壓力云圖Fig.7 Pressure contour of cruise shape

表1 不同展向位置的機(jī)翼相對(duì)扭角Table 1 Relative torsional angles of the wing at different spanwise sections

利用文中建立的相關(guān)技術(shù)對(duì)型架外形進(jìn)行靜氣彈計(jì)算，本文主要驗(yàn)證了型架構(gòu)型在巡航狀態(tài)下氣動(dòng)特性。在利用稀疏表面網(wǎng)格節(jié)點(diǎn)與空間網(wǎng)格頂點(diǎn)構(gòu)建Delaunay圖時(shí)，Delaunay單元數(shù)為6090個(gè)，計(jì)算時(shí)間為30s，之后由稀疏表面網(wǎng)格節(jié)點(diǎn)與空間網(wǎng)格頂點(diǎn)進(jìn)行映射單元變形。本文采用以上方法既避免了Delaunay單元構(gòu)建輸入點(diǎn)的復(fù)雜操作，又從很大程度上降低了RBF大型矩陣的求解量，從而大幅度提高網(wǎng)格變形效率。表2給出了相同網(wǎng)格條件下（以DLR－F6翼身組合體為對(duì)象），不同變形網(wǎng)格技術(shù)與本文方法工作效率的對(duì)比，可以看出本文提出的新型變形網(wǎng)格方法具備較高的網(wǎng)格更新速度。圖10給出了經(jīng)過靜氣動(dòng)彈性計(jì)算后不同網(wǎng)格區(qū)域頂點(diǎn)位置在RBF徑向基函數(shù)技術(shù)操作下的更新。圖11進(jìn)一步給出了空間頂點(diǎn)更新后對(duì)應(yīng)Delaunay圖空間截面。在該計(jì)算模型靜氣彈計(jì)算過程中，RBF徑向基函數(shù)技術(shù)對(duì)330個(gè)不同網(wǎng)格區(qū)域1000個(gè)頂點(diǎn)的操作僅需要0.05s時(shí)間，而對(duì)于頂點(diǎn)更新后Delaunay圖單元變形進(jìn)行的空間1300萬網(wǎng)格節(jié)點(diǎn)變形僅需要1.5s，可見本文構(gòu)建的徑向基函數(shù)與Delaunay圖映射網(wǎng)格變形技術(shù)效率較高，比較適用于大型工程靜氣動(dòng)彈性計(jì)算中。整個(gè)靜氣動(dòng)彈性經(jīng)過10次迭代基本收斂。圖12給出了型架構(gòu)型巡航外形與巡航設(shè)計(jì)外形的對(duì)比，可以看出兩者基本完全重合，驗(yàn)證了修正方法的正確性與靜氣動(dòng)彈性分析的可靠性與較高效率。表3給出了兩者在氣動(dòng)特性方面的對(duì)比，兩者氣動(dòng)特性接近，誤差較小。圖13～圖16給出了展向典型站位壓力分布對(duì)比，壓力分布形態(tài)相同基本重合。

表2 不同變形網(wǎng)格技術(shù)效率對(duì)比Table 2 Efficient comparison between different grid deformation methods

圖8 型架外形與巡航外形對(duì)比Fig.8 Jig shape and cruise shape

圖9 型架外形與巡航外形局部視圖Fig.9 Local view of jig shape and cruise shape

圖10 RBF操作下區(qū)域頂點(diǎn)分布Fig.10 Block vertexs distribution under RBF operation

圖11 Delaunay四面體空間截面Fig.11 Delaunay space section

圖12 巡航設(shè)計(jì)外形與型架巡航外形對(duì)比Fig.12 Comparision of cruise designed shape and jig cruise shape

表3 巡航設(shè)計(jì)外形與型架巡航外形氣動(dòng)特性對(duì)比Table 3 Comparison of aerodynamic performance between cruise design shape and jig shape at cruise condition

圖13 y／b＝0.19展向位置壓力分布Fig.13 Pressure distribution of y／b＝0.19spanwise section

圖14 y／b＝0.3展向位置壓力分布Fig.14 Pressure distribution of y／b＝0.3spanwise section

圖15 y／b＝0.5展向位置壓力分布Fig.15 Pressure distribution of y／b＝0.5spanwise section

圖16 y／b＝0.8展向位置壓力分布Fig.16 Pressure distribution of y／b＝0.8spanwise section

5 結(jié)束語

本文研究了采用RBF徑向基函數(shù)、Delaunay圖映射技術(shù)的網(wǎng)格變形方法。研究建立了基于靜氣動(dòng)彈性計(jì)算的飛行器型架外形修正設(shè)計(jì)、驗(yàn)證方法。對(duì)某型支線客機(jī)進(jìn)行了型架外形逆推、修正及靜氣動(dòng)彈性驗(yàn)證，得出如下結(jié)論：

（1）基于RBF徑向基函數(shù)、Delaunay圖映射的網(wǎng)格變形技術(shù)充分利用多塊網(wǎng)格頂點(diǎn)，在避免Delaunay圖單元構(gòu)建時(shí)輸入點(diǎn)的復(fù)雜操作的同時(shí)充分利用RBF對(duì)網(wǎng)格頂點(diǎn)的高效操作，與傳統(tǒng)變形網(wǎng)格技術(shù)相比，計(jì)算效率較高，適應(yīng)大幅度網(wǎng)格變形，適合于大型工程靜氣動(dòng)彈性計(jì)算。

（2）逆推的飛行器型架構(gòu)型表面光滑，變形合理，說明所采用“映射曲面”CFD／CSD數(shù)據(jù)交換技術(shù)插值精度高。

（3）通過對(duì)型架外形的氣動(dòng)彈性計(jì)算，驗(yàn)證了文中型架構(gòu)型修正的正確性以及靜氣動(dòng)彈性分析方法的可靠性、高效性。

［1］XIE C C，WU Z G，YANG C.Aeroelastic analysis of flexible large aspect ratio wing［J］.Journal of Beijing University of Aeronautics and Astronautics，2003，29（12）：1087－1090.（in Chinese）謝長川，吳志剛，楊超.大展弦比柔性機(jī)翼的氣動(dòng)彈性分析［J］.北京航空航天大學(xué)學(xué)報(bào)，2003，29（12）：1087－1090.

［2］GUAN D.Aeroelastic experiment［M］.Beijing Aeronautics College press，1986.（in Chinese）管德.氣動(dòng)彈性試驗(yàn)［M］.北京航空學(xué)院出版社，1986.

［3］YANG G W，ZHENG G N.Aircraft jig shape correction method based on static aeroelastic analyses［J］.Advances in Aeronautical Science and Engineering，2011，2（2）：113－150.（in Chinese）楊國偉，等.基于靜氣動(dòng)彈性效應(yīng)的飛機(jī)型架外形修正方法研究［J］.航空工程進(jìn)展，2011，2（2）：113－150.

［4］ZHAN H，CHENG S X，ZHU J.An effective aerodynamic and aeroelastic coupled design method for complicated wing shape［J］.Journal of Northwestern Polytechnical University，2009，27（1）：100－104.（in Chinese）詹浩，等.考慮氣動(dòng)彈性影響的機(jī)翼復(fù)雜氣動(dòng)外形設(shè)計(jì)研究［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2009，27（1）：100－104.

［5］ALLEN C B.Unified approach to CFD－CSD interpolation and mesh motion using radial basis functions［R］.AIAA Paper 2007－3084，2007.

［6］LIU X Q，QIN N.Fast dynamic grid deformation based on Delaunay graph mapping［J］.Journal of Computational Physics，2006，211：405－423.

［7］BUHMANN M.Radial basis functions［M］.Cambridge University Press，2005.

［8］CIZMAS P，GARGOLOFF J.Mesh generation and deformation algorithm for aeroelasticity simulations［R］.AIAA 2007－556.45th Aerospace Sciences Meeting［C］.Reno，NV，2007.

［9］PARK J，SANDBERG I W.Universal approximation using radialbasis－function networks［J］.Neural Computation，1991.

［10］LIGHT W A.Some aspects of radial basis function approximation.approximation theory，spline functions and applications［J］.NATO ASI Series，1992，356：163－190.

［11］KARIM A，ADELI H.Radial basis function neural network for work zone capacity and queue estimation［J］.J.Transp.Eng.，2003，129（5）：494－503.

［12］ZHANG X，SONG K Z，LU M W，et al.Meshless methods based on collocation with radial basis functions［J］.Computational Mechanics，26（4）：333－343.

［13］SARLER B.A radial basis function collocation approach in computational fluid dynamics［J］.CMES，2005，7（2）：185－193.

［14］WENDLAND H.Fast evaluation of radial basis functions：Methods based on partition of Unity［M］／／Approximation Theory X：Wavelets，Splines，and Applications.Nashville，Texas，USA：Vanderbilt University Press，2002：473－483.

［15］LI L Z.Turbine blade temperature transfer using the load surface method［J］.Computer Aided Design，2007，39：494－505.

［16］ZHU X X.Free curves and surfaces modeling technology［M］.Science Press，2000.（in Chinese）朱心雄.自由曲線曲面造型技術(shù)［M］.科學(xué)出版社，2000.

［17］ZHAO Y H.Aeroelastic dynamics and control［M］.Science Press，2007.（in Chinese）趙永輝.氣動(dòng)彈性力學(xué)與控制［M］.科學(xué)出版社，2007.

［18］SHEN K Y，GUAN D.The principle of elastic mechanics［M］.Shanghai science and Technology Press，1982.（in Chinese）沈克揚(yáng)，管德.彈性力學(xué)原理［M］.上?？茖W(xué)技術(shù)文獻(xiàn)出版社，1982.