樣條有限元條法

1、第6章 樣條有限元條法§6.1 緒論20世紀(jì)60年代在R.W.Clough命名后迅速發(fā)展起來(lái)的有限單元法，被公認(rèn)為結(jié)構(gòu)分析最強(qiáng)有力的工具，理論上可適用于所有的結(jié)構(gòu)。但對(duì)于規(guī)則區(qū)域的結(jié)構(gòu)，有限元法的解題效率不如差分法；在計(jì)算機(jī)技術(shù)取得突破性進(jìn)展之前，解決大型空間結(jié)構(gòu)工程設(shè)計(jì)問(wèn)題，往往需要大型計(jì)算機(jī)和昂貴的費(fèi)用，使有限元法受到種種限制；有些結(jié)構(gòu)問(wèn)題用有限元法也難以解決，需要?jiǎng)?chuàng)造一些新的數(shù)值方法，于是又產(chǎn)生了有限條法、邊界元法、樣條元法以及新的加權(quán)殘數(shù)（余量）法等等。適用于幾何形狀規(guī)則和邊界條件簡(jiǎn)單的土建結(jié)構(gòu)的一種行之有效的特殊有限元法半解析有限條法，是由Y.K.Cheung（張佑啟）和E

2、.L.Wilson在1968年提出的17。與有限元法相比，主要不同點(diǎn)在于所取得位移函數(shù)，一般是以多項(xiàng)式和正交級(jí)數(shù)乘積的形式給出，使得彈性力學(xué)問(wèn)題降維，從而使總剛度矩陣大大降階，既省機(jī)時(shí)，精度也高。但對(duì)于非簡(jiǎn)支端邊界條件，級(jí)數(shù)耦合，計(jì)算繁雜，尤其是對(duì)于集中荷載作用和內(nèi)部支承情況，所取項(xiàng)數(shù)多，收斂慢。對(duì)于處理沿跨向材料突變或是變截面問(wèn)題，結(jié)果也不能令人滿意。1979年石鐘慈提出樣條有限元法18，用三次B樣條變分方法解規(guī)則區(qū)域上板梁組合彈性結(jié)構(gòu)的平衡問(wèn)題，導(dǎo)出適用于各種邊界條件的統(tǒng)一計(jì)算格式，比通常的有限元法計(jì)算量少，精度高，便于在小型計(jì)算機(jī)上實(shí)現(xiàn)。稍后，秦榮提出以樣條函數(shù)、梁振動(dòng)函數(shù)及能量變分為基

3、礎(chǔ)的樣條有限點(diǎn)法和后來(lái)的樣條子域法19。1982年Y.K.Cheung等人又提出結(jié)構(gòu)分析的樣條有限條法，以克服經(jīng)典有限條法的缺點(diǎn)。作者在分析空間懸掛結(jié)構(gòu)吊橋時(shí)，采用元條相結(jié)合的方法20，將主索剛度疊加于橋面板條單元的結(jié)線上，分析索結(jié)構(gòu)的動(dòng)力特性取得成果。§6.2 B樣條與樣條函數(shù)樣條函數(shù)是現(xiàn)代函數(shù)逼近的一個(gè)十分活躍的分支，是計(jì)算方法的一個(gè)重要基礎(chǔ)，已得到廣泛應(yīng)用。樣條函數(shù)來(lái)源于實(shí)際生產(chǎn)中的樣條曲線，由于樣條函數(shù)是一個(gè)分段多項(xiàng)式，利用它去逼近任意函數(shù)，具有更大的靈活性和適應(yīng)性。樣條函數(shù)可以利用基本樣條（如B樣條）函數(shù)的線性組合來(lái)構(gòu)造，在計(jì)算力學(xué)中，常用B樣條函數(shù)來(lái)構(gòu)造位移函數(shù)及應(yīng)力函數(shù)

4、。 B樣條及其性質(zhì)B樣條函數(shù)是基本樣條函數(shù)（簡(jiǎn)稱B樣條），它與廣義函數(shù)函數(shù)有密切的內(nèi)在聯(lián)系，并可以逼近函數(shù)，所以也稱樣條函數(shù)。B樣條函數(shù)系可作成樣條函數(shù)空間的一組基。n次B樣條函數(shù)可寫成統(tǒng)一的表達(dá)式 (）式中，是二項(xiàng)式系數(shù)的組合表示，樣條結(jié)點(diǎn) 。吊橋動(dòng)力分析主要應(yīng)用一次和三次B樣條，其分段表達(dá)式和曲線形狀如下： (） () (a) 一次B樣條 (b) 三次B樣條圖 B樣條曲線n次B樣條函數(shù)具有緊致性（即n(x)及其導(dǎo)數(shù)以區(qū)間為緊致支撐集，此外，處處為零）和分段光滑性（為n次分段多項(xiàng)式，在（-，）區(qū)間，具有n-1 階連續(xù)導(dǎo)數(shù)）以及對(duì)稱（即n(-x)=n(x)）和面積不變（亦即）等特性。尤其是奇次

5、樣條函數(shù)具有極小模性質(zhì)（可用來(lái)量度變形能極小性質(zhì)）以及最佳逼近性質(zhì)，顯示了很大的數(shù)學(xué)與力學(xué)的使用價(jià)值。 樣條函數(shù)樣條函數(shù)曲線來(lái)源于生產(chǎn)實(shí)踐中的數(shù)學(xué)放樣。定義一次樣條函數(shù)為二階廣義微分方程 的解，即 (）式中，xj為結(jié)點(diǎn)，下標(biāo)“+”為截?cái)喾?hào)，即定義： (）可見(jiàn)，一次樣條函數(shù)是連續(xù)，且分段線性的折線函數(shù)，其力學(xué)意義是集中力作用下的彈性弦，它表征了纜索受吊桿作用力時(shí)的實(shí)際情況。我們也可以定義三次樣條函數(shù)為四階廣義微分方程 的解，其一般表達(dá)式為 ()顯然，這是一個(gè)分段三次多項(xiàng)式（當(dāng)j=0時(shí)，便退化為普通的三次多項(xiàng)式），有二階連續(xù)導(dǎo)數(shù)，三階階躍導(dǎo)數(shù)。如果與梁的撓曲線方程y(4)(x)=q(x)對(duì)照，可

6、知三次樣條函數(shù)相當(dāng)于彈性梁在結(jié)點(diǎn)xj處受集中荷載j作用下的撓曲線。因此，可用以反映吊橋加勁梁受吊桿力作用的實(shí)際情況。同樣，我們可以定義更高次的樣條函數(shù)，例如，在研究連續(xù)應(yīng)力場(chǎng)和應(yīng)變場(chǎng)時(shí)，定義了五次樣條函數(shù)21。二次樣條函數(shù)也同樣有其力學(xué)背景。這里不贅述。根據(jù)展開(kāi)定理，n次樣條函數(shù)Sn(x)可表為n次B樣條函數(shù)n(x)及其平移n(x-k)的線性組合，即 (）當(dāng)n為奇數(shù)時(shí)，若在0,的區(qū)間上取N等分劃，則h=/N，在第k結(jié)點(diǎn)，有xk=kh，于是上式可表為 ()對(duì)于一次和三次B樣條，有 (）和 (）而 上以為結(jié)點(diǎn)的一次B樣條的一組基函數(shù)，張成N+1維空間。類似的，為三次B樣條的一組基函數(shù)，張成N+3維

7、空間。在線性空間等價(jià)的前提下，可改換基函數(shù)的形式，構(gòu)造出適合結(jié)構(gòu)位移邊界條件的基函數(shù)，使統(tǒng)一處理各種邊界條件得到方便。例如，將（）式改寫后表示位移函數(shù) (）由于基函數(shù)的局部緊湊，若要計(jì)算某樣條結(jié)點(diǎn)xi的位移值，上式中最多只有三項(xiàng)不為零，即 (）在求出廣義參數(shù)后，經(jīng)上式計(jì)算，才得到待求的位移值（與有限元不同）。由于k(x)在結(jié)點(diǎn)上的值是一些簡(jiǎn)單的現(xiàn)成的數(shù)，故需要的計(jì)算量很小。這些“現(xiàn)成數(shù)”已分幾種邊界情況編入子程序中。是一個(gè)與三次B樣條有關(guān)的基函數(shù)，考慮常見(jiàn)的邊界條件，寫成： 以上的基函數(shù)對(duì)應(yīng)于自由邊界。若為簡(jiǎn)支邊界，則刪去第一行（列）或末一行（列）；若為固定邊界，則刪去首兩行（列）或末兩行（列

8、）。若結(jié)構(gòu)對(duì)稱，可簡(jiǎn)化取半計(jì)算，將半結(jié)構(gòu)分成N等分劃，根據(jù)滑動(dòng)固支邊界（W=0），僅刪去對(duì)稱軸一端的三行三列，代以新的二行二列，即 反對(duì)稱時(shí)，取半結(jié)構(gòu)的對(duì)稱軸一端仍按鉸支處理。 利用參數(shù)法19，還可構(gòu)造出與三次B樣條函數(shù)有關(guān)的新的基函數(shù)，得到具有物理意義的位移參數(shù)和廣義參數(shù)混合的位移函數(shù)，使處理位移邊界條件和連續(xù)條件都很方便。經(jīng)過(guò)對(duì)B樣條函數(shù)的求導(dǎo)和積分，計(jì)算出系數(shù)矩陣提供編程使用。 B樣條函數(shù)的求導(dǎo)和積分B樣條函數(shù)的求導(dǎo)和積分公式為 ()式中，當(dāng)j為正整數(shù)時(shí)，表示求j次導(dǎo)數(shù)；當(dāng)j為負(fù)整數(shù)時(shí)，表示積分(-j)次；當(dāng)j=0時(shí)，即為（）式。以B樣條基函數(shù)（或其導(dǎo)數(shù)）為元素的矩陣乘積的積分是一組很有

9、實(shí)用意義的系數(shù)矩陣，其具體形式可事先算出，列表待查。吊橋的動(dòng)力分析將用到以下的系數(shù)矩陣，記為： ()以I4y為例（暫不考慮修改邊界），說(shuō)明矩陣乘積的積分過(guò)程。 令 則 ，反復(fù)進(jìn)行分部積分。由于 于是 ()同理，可進(jìn)行其余矩陣的積分。特殊的系數(shù)矩陣（如修改邊界時(shí)，可能是3的某種線性組合）根據(jù)問(wèn)題的需要區(qū)別處理，隨時(shí)算出。由于三次B樣條的緊湊性，所生成的系數(shù)矩陣一般是半帶寬為4的帶狀矩陣。§6.3 樣條有限元條的位移模式 吊橋簡(jiǎn)化假定以吊橋?yàn)槔诶脴訔l有限元條法進(jìn)行吊橋動(dòng)力特性分析時(shí)，作如下假定： 小振幅振動(dòng)，所有材料遵循虎克定律。 纜索忽略彎曲剛度。 初始靜載全部由主索承受，索呈二

10、次拋物線的理想索形。 忽略細(xì)長(zhǎng)加勁橋面結(jié)構(gòu)的剪切變形和轉(zhuǎn)動(dòng)慣量影響。 主索的重量沿跨長(zhǎng)均布，各跨加勁結(jié)構(gòu)等截面。 吊桿不可伸長(zhǎng)，在豎直扭轉(zhuǎn)振動(dòng)時(shí)保持垂直。 振動(dòng)引起的索的附加水平張力和塔頂?shù)母郊訅毫珊雎圆挥?jì)。 忽略加勁結(jié)構(gòu)的初應(yīng)力和初始曲率。 結(jié)構(gòu)離散與位移模式根據(jù)上述的線性化假定，僅考慮橋面結(jié)構(gòu)的彎扭振型時(shí)，懸掛的空間結(jié)構(gòu)被離散為索單元，塔單元和橋面殼條元體系。1. 索作為一維樣條弦桿單元考慮，采用一次B樣條（由于實(shí)際吊桿密布，有時(shí)索也采用三次B樣條，以便統(tǒng)一處理）。僅在豎直平面內(nèi)有一個(gè)自由度（WC）。由剛性吊桿的假定，其有Wci=Wi，這里Wi為加勁橋面的豎向位移。索與橋面系在順橋向都取

11、相同的N個(gè)等分劃，即有N+1個(gè)結(jié)點(diǎn)。位移模式為 ，而，矩陣中各元素i(y)是與B樣條函數(shù)有關(guān)的基函數(shù)。Wci=W-1 W0 W1WiWN+1 T ,其中Wci是待定的與時(shí)間有關(guān)的豎向結(jié)點(diǎn)位移參數(shù)。2. 橋面系視為一般的正交異性折板結(jié)構(gòu)（橋?qū)挷淮髸r(shí)，也可視為一維梁），包括了板梁、箱梁和閉合桁架等典型橋面結(jié)構(gòu)，并進(jìn)一步離散為樣條有限條的矩形殼條單元（如圖所示），坐標(biāo)軸x,y為其主方向。順橋向采用N個(gè)等分劃的三次B樣條，橫橋向采用二結(jié)點(diǎn)低階條的常規(guī)形函數(shù)17。假定不計(jì)殼條彎曲和薄膜效應(yīng)之間的相互作用，每條結(jié)線（內(nèi)設(shè)結(jié)點(diǎn)）上有四個(gè)自由度，用位移參數(shù)u,v,w,表示。特殊情況作為板條單元，只計(jì)彎曲效應(yīng)，

12、每條結(jié)線只有兩個(gè)自由度，位移參數(shù)為w,(為埃米特插值所需)。位移模式表為： ()式中，A是形函數(shù)矩陣，而是位移參數(shù)列矩陣。其中 (）矩陣元素Ni(x)為低階條的常規(guī)形函數(shù)，其中： 子矩陣 式中，加“”記號(hào)的基函數(shù)一般是根據(jù)邊界條件局部修改過(guò)的。 子矩陣 圖 樣條矩形殼條單元 3. 橋塔被離散為四個(gè)或更多的在豎直面內(nèi)撓曲的懸臂梁（柱）單元，采用三次B樣條，取N個(gè)等分劃，每個(gè)內(nèi)結(jié)點(diǎn)僅由一個(gè)自由度（Vt）,位移模式為：，這里，和Vtj的定義類似于索單元。§6.4 動(dòng)力特性矩陣樣條有限元條法基本上是變分問(wèn)題的直接法（Ritz-Galerkin法）和以位移逼近的有限元法的綜合。 結(jié)構(gòu)振動(dòng)的泛函

13、 根據(jù)哈密頓原理，可知結(jié)構(gòu)振動(dòng)的泛函為 ()式中：t為時(shí)間，t1tt2；U為結(jié)構(gòu)的應(yīng)變能；V為結(jié)構(gòu)的外力功；T為結(jié)構(gòu)的動(dòng)能。設(shè)結(jié)構(gòu)振動(dòng)時(shí)的位移函數(shù)為ft =u，v，wT,為板殼單位面積質(zhì)量，t為應(yīng)變，t為曲率（扭率），qt為分布外動(dòng)荷載，則 () () ()研究無(wú)阻尼自由振動(dòng)時(shí)，動(dòng)荷載qt=0。由于微幅諧振動(dòng)，振動(dòng)的位移函數(shù)是一個(gè)周期函數(shù)，在某一周期T=2/內(nèi)對(duì)時(shí)間t積分，可得結(jié)構(gòu)自由振動(dòng)的泛函為 ()代入吊橋結(jié)構(gòu)的位移模式 ()令 ()把（）式寫成，由變分原理轉(zhuǎn)化為函數(shù)的極值問(wèn)題，即 ，得到吊橋不計(jì)阻尼自由振動(dòng)的動(dòng)力矩陣方程： ()求解該廣義特征問(wèn)題，得出結(jié)構(gòu)的固有頻率及其相應(yīng)的振型。 動(dòng)力

14、特性矩陣 在給出殼條或單元在平衡位置的勢(shì)能和動(dòng)能的表達(dá)式后，由影響系數(shù)的定義，容易得到條和元的剛度矩陣s和一致質(zhì)量矩陣m。索的剛度、質(zhì)量與條的剛度、質(zhì)量在結(jié)線上疊加；塔的剛度、質(zhì)量可作為獨(dú)立的矩陣塊。通過(guò)坐標(biāo)變換和根據(jù)邊界條件局部修改樣條函數(shù)后，集成結(jié)構(gòu)總剛度矩陣K和總一致質(zhì)量矩陣M。對(duì)于矩形扁殼條，廣義應(yīng)變?yōu)槠渲校瑧?yīng)變矩陣： 廣義應(yīng)力為 對(duì)于正交異性材料，彈性矩陣： 其中: 以上各式中，Ex、Ey、x、y和G是彈性常數(shù)，t是條的厚度。當(dāng)各向同性時(shí)，Ex=Ey=E，。殼條應(yīng)變能： （1）殼條剛度矩陣 ()上式各分塊矩陣中，s12=s21= 0 其中：以上各式中的系數(shù)矩陣，僅I5y為反對(duì)稱矩陣，

15、其余均為對(duì)稱矩陣。（2）殼條的一致質(zhì)量矩陣殼條動(dòng)能：式中，是單位體積的重量，t為板厚，A=N。故殼條的一致質(zhì)量矩陣為 ()各分塊矩陣為： m12=m21= 0， ， 橋面系沿橫向可以有不同方位布置的條元，如圖所示，條元的局部坐標(biāo)系y軸和結(jié)構(gòu)整體坐標(biāo)系y軸平行（相鄰條交線重合），為x和x軸之間的夾角（順時(shí)針為正）。殼條坐標(biāo)轉(zhuǎn)換矩陣為 (6.4.11)其中 圖 坐標(biāo)轉(zhuǎn)換 （3）索的彈性剛度矩陣將索的位移模式代入索張力變化所儲(chǔ)存應(yīng)變能表達(dá)式，導(dǎo)出索的彈性剛度矩陣為 ()式中，為索的等效長(zhǎng)度，為傾角，為一系數(shù)矩陣，為單位長(zhǎng)度的索重加橋面結(jié)構(gòu)重，為索的靜水平張力。 （4）索的重力剛度矩陣索的動(dòng)伸長(zhǎng)對(duì)索撓

16、度的影響與索豎向振幅相比是高階微量，可忽略不計(jì)。代入索的位移模式和相容方程，索的重力勢(shì)能可表達(dá)為 則索的重力剛度矩陣為 () （5）索的一致質(zhì)量矩陣索的動(dòng)能表達(dá)式： 則索的一致質(zhì)量矩陣為 () （6）塔柱的剛度矩陣塔柱取等效系統(tǒng)分析，假設(shè)其彎曲和扭轉(zhuǎn)由各柱在豎直面內(nèi)的撓曲形成?？紤]塔柱的應(yīng)變能和外力勢(shì)能，得到塔柱彈性剛度矩陣為 ()式中，I1Z（及下式I3Z、I4Z）為一系數(shù)矩陣，參照I1Y等。塔柱變剛度時(shí)，若慣性矩It為z的連續(xù)函數(shù)，則可參與積分，給出新的系數(shù)矩陣；若剛度呈階梯形變化，則應(yīng)分段，各等分劃處理。塔柱的幾何剛度矩陣考慮塔頂?shù)撵o、動(dòng)壓力PW和P(t)以及塔頂?shù)刃椈蓜偠萲,寫成 (

17、) （7）塔柱的一致質(zhì)量矩陣考慮單塔柱的動(dòng)能，得到單塔柱的一致質(zhì)量矩陣表達(dá)式： ()對(duì)于三跨簡(jiǎn)支吊橋，如果橋面結(jié)構(gòu)的材料、構(gòu)造都相同，可按跨分段,各取等分劃。根據(jù)簡(jiǎn)支邊界條件局部修改B樣條后，三跨的基函數(shù)集成對(duì)角分塊矩陣，其階數(shù)可隨各跨分劃數(shù)的不同而不同，樣條矩形殼條的剛度矩陣形式同單跨，系數(shù)矩陣按跨分塊，集成對(duì)角帶狀矩陣；索的彈性剛度矩陣也由三跨分塊集成，對(duì)稱局部滿陣。位移參數(shù)列陣也由各跨集成。對(duì)于三跨連續(xù)的吊橋橋面結(jié)構(gòu)，可在中間支承處分段。為方便起見(jiàn)，系數(shù)矩陣需經(jīng)修改，構(gòu)造出具有混合參數(shù)的位移函數(shù)，在此不贅述。 §6.5 算例 例1 用樣條有限元（按一維工程梁）計(jì)算如圖所示單跨吊

18、橋的豎向彎曲振動(dòng)。特性數(shù)據(jù)如下：加勁梁：L=853.44m, W=4246.016/m, ES=203932.16MN/m2, IS=7.695953m4 索：f=70.714m, AC=0.123548m2, EC=179129.6MN/m2, LE=1219.2m, HW=53578kN.圖 單跨吊橋?qū)τ谝痪S梁?jiǎn)卧O(shè)W=w剛度矩陣: 一致質(zhì)量矩陣：幾種方法計(jì)算結(jié)果示于表。比較表明，樣條有限元法具有較高的計(jì)算精度。表 單跨吊橋豎向振動(dòng)頻率頻率序自 振 頻 率 （rad/sec）對(duì) 稱 振 型反 對(duì) 稱 振 型解析解有限元S-12S-16解析解有限元S-12S-1611.4001.39751

19、.39741.39741.33181.33301.33041.330422.6962.70472.70222.70224.49014.48704.48854.484936.84726.85656.84369.71399.71639.76289.7073413.118713.270413.136217.020517.046217.456017.0376521.515322.456321.683826.412826.563628.360126.5569注：表中S-12和S-16分別表示樣條有限元條法取12分劃和16分劃。例2 三跨簡(jiǎn)支吊橋舊Tacoma橋（該橋于1940年在不到20m/s的風(fēng)速下發(fā)生顫振破壞）主要設(shè)計(jì)數(shù)據(jù)如下：ES=2.038318