彈性壓桿穩(wěn)定問題的ODE解法

徐杏華，李 朝

(1.孝感學(xué)院城市建設(shè)學(xué)院，湖北 孝感 432000，2.長春市軌道交通集團(tuán)有限公司，吉林長春 130012)

在結(jié)構(gòu)設(shè)計(jì)中，需要對結(jié)構(gòu)進(jìn)行強(qiáng)度驗(yàn)算和穩(wěn)定驗(yàn)算，然而對于以受壓為主的結(jié)構(gòu)來說，其結(jié)構(gòu)更容易失穩(wěn)。壓桿作為結(jié)構(gòu)失穩(wěn)的典型代表，它的失穩(wěn)，輕則引起構(gòu)件失效，重則引起整個(gè)結(jié)構(gòu)的破壞，造成嚴(yán)重的事故。在壓桿穩(wěn)定問題研究方面，楚中毅等［1，2］對穩(wěn)定問題的精確解法做了一些探討，龍馭球［3］也從不同角度總結(jié)了穩(wěn)定問題的各種解法。然而對于無限自由度桿系結(jié)構(gòu)線性穩(wěn)定問題，解析法求解的精度往往與其計(jì)算時(shí)所選取的撓曲線的函數(shù)形式有很大的關(guān)系，利用普通結(jié)構(gòu)力學(xué)的方法根本不能得到精確解［4］，而采用常規(guī)有限元法，往往需要對網(wǎng)格進(jìn)行細(xì)分，這樣做不僅計(jì)算量大、計(jì)算效率低，而且結(jié)果也不夠精確。

隨著常微分方程(即Ordinary Differential Equation，簡記作ODE)數(shù)值解法的發(fā)展，尤其是近10年來一系列常微分方程求解器(即 Ordinary Differential Equation Solver)通用程序相繼問世［5，6，7］，使直接針對結(jié)構(gòu)穩(wěn)定問題的數(shù)值解析法成為可能。對于一般的平面桿系結(jié)構(gòu)，常微分方程求解器不僅可以給出精確的臨界荷載和相應(yīng)的失穩(wěn)變形形態(tài)，而且還可以給出高階失穩(wěn)荷載和形態(tài)。本文推導(dǎo)出了無限自由度彈性壓桿穩(wěn)定問題的控制微分方程，算例結(jié)果與解析解以及常規(guī)有限元解的比較表明該方法的求解精度和效率較高。

1 壓桿的撓曲線方程

1.1 控制微分方程的推導(dǎo)

考慮圖1所示的最具一般性的壓桿，其支撐端可以同時(shí)受三個(gè)力作用，即:壓力P，約束反力R和約束力偶矩M0，相應(yīng)的彎矩方程應(yīng)為:

圖1 壓桿受力狀況Fig.1 Stress state of compression bar

1.2 控制方程的標(biāo)準(zhǔn)化

由于常微分方程求解器是按照標(biāo)準(zhǔn)的ODE形式研制的，只具有求解標(biāo)準(zhǔn)的非線性常微分方程問題的功能，不具有直接求解特征值問題的功能［5］。因此對于上述關(guān)于邊值問題的常微分方程組的特征值問題，應(yīng)該利用一些ODE變換技巧將其轉(zhuǎn)換為COLSYS［7，8］所能接受的標(biāo)準(zhǔn)的非線性常微分方程形式。為此需做以下工作:

區(qū)間映射:利用區(qū)間映射技巧進(jìn)行坐標(biāo)變換，將非規(guī)則的求解區(qū)間映射到確定的標(biāo)準(zhǔn)單位區(qū)間，即將待定的特征值問題定義在標(biāo)準(zhǔn)區(qū)間［0，1］上，使結(jié)構(gòu)的每一個(gè)單元結(jié)點(diǎn)處的坐標(biāo)都為0或1(如圖2)。例如:對于一長度為L的單元，我們可令

圖2 坐標(biāo)轉(zhuǎn)換圖Fig.2 Diagram of coordinate transformation

經(jīng)坐標(biāo)變換后，此變截面壓桿的撓曲線方程可轉(zhuǎn)化為:

式中λ=P，Y=Y(ξ)，()＇=d()/dξ，n表示因結(jié)構(gòu)剛度或質(zhì)量變化而將結(jié)構(gòu)劃分的段數(shù)。

(2)構(gòu)造平凡ODE:為將特征值轉(zhuǎn)化為標(biāo)準(zhǔn)的非線性O(shè)DE問題，需要為待定的特征值建立一個(gè)平凡ODE，從而保證λ為常數(shù)的前提下，將它引到ODE中求解。即

(3)構(gòu)造等價(jià)ODE:將一重積分問題轉(zhuǎn)換為一階常微分方程問題，以便于COLSYS求解，即取振型歸一化條件為

這一歸一化條件相當(dāng)于將Rayleigh商中的分母取為單位值，對上式(4)進(jìn)行坐標(biāo)變換后，有

再利用等價(jià)的ODE技巧將以上積分式轉(zhuǎn)換為標(biāo)準(zhǔn)的ODE問題

這樣既可以計(jì)算定積分求解上式，又可以將某些積分條件(如歸一化條件)化為等價(jià)的ODE提法而加入原問題求解。于是式(3)、式(4)、式(7)便形成了一組標(biāo)準(zhǔn)的非線性的常微分方程組，在進(jìn)行計(jì)算時(shí)，只要用戶提供一個(gè)適當(dāng)?shù)某跏冀猓涂芍苯永脴?biāo)準(zhǔn)的常微分方程求解器的非線性功能進(jìn)行求解。

2 求解器解法程序結(jié)構(gòu)及其調(diào)用

常微分方程求解器主要是通過調(diào)用常微分方程邊值問題通用程序COLSYS來實(shí)現(xiàn)其求解功能的，它是目前最為可靠的常微分方程邊值問題求解器之一，也是面向常微分方程常用的支撐軟件。COLSYS是以子程序的形式由用戶進(jìn)行調(diào)用，在調(diào)用時(shí)只需要用戶為COLSYS提供某些參數(shù)。例如:每個(gè)子區(qū)間的配置點(diǎn)數(shù)，初始區(qū)間網(wǎng)格劃分的子區(qū)間數(shù)，給定誤差限的解分量、誤差和誤差限值以及導(dǎo)數(shù)的數(shù)目，初始網(wǎng)格的劃分，非限性問題初始解的提供形式，非線性問題是否是敏感型等。

控制微分方程及其相應(yīng)的邊界與連接條件的編程求解是工作的核心部分，這個(gè)模塊需要用戶根據(jù)具體的問題在外部定義FSUB、DFSUB、GSUB、DGSUB、SOLUTN五個(gè)子程序以供調(diào)用。下面將詳細(xì)地介紹這五個(gè)子程序的功能:

(1)FSUB是一個(gè)提供常微分方程信息的子程序名。值得注意的是在進(jìn)行編程求解前，應(yīng)先利用ODE變換技術(shù)將這些常微分方程轉(zhuǎn)換為能夠?yàn)镃OLSYS所接受標(biāo)準(zhǔn)的非線性常微分方程形式，以便使結(jié)構(gòu)的每一個(gè)單元結(jié)點(diǎn)處的坐標(biāo)都為0或1。

(2)DFSUB是一個(gè)計(jì)算FSUB中常微分方程的雅可比矩陣的子程序名。

(3)GSUB是一個(gè)提供邊界條件與連接條件的子程序名。

(4)DGSUB是一個(gè)計(jì)算GSUB中Gi(X，Z)的雅可比矩陣的子程序名。

(5)SOLUTN是一個(gè)計(jì)算非線性問題初始解的子程序名。只有求解非線性問題時(shí)才會(huì)用到此子程序。

除以上5個(gè)外部子程序外，還需編制一些輔助子程序，以及主程序和子程序之間的聯(lián)系程序單元即接口塊，這樣一個(gè)完整的調(diào)用COLSYS求解器解決ODE問題的程序結(jié)構(gòu)即告完成，可以進(jìn)行線性或非線性O(shè)DE體系問題的編程計(jì)算。

3 算例與計(jì)算結(jié)果分析

計(jì)算圖3所示為一懸臂壓桿的臨界荷載和失穩(wěn)模態(tài)。其中彈性模量E=1，截面慣性矩I=1，桿長l=1，初始軸力P=10。該問題的理論解答為Pcr=π2EI/(2l)2=2.46740110，表1為本文方法的計(jì)算結(jié)果，表2為常規(guī)有限元法的計(jì)算結(jié)果。

圖3 懸臂壓桿及其失穩(wěn)模態(tài)Fig.3 Cantilever bar and corresponding instability mode

表1 ODE解計(jì)算結(jié)果Table 1 Result of ODE method

表2 常規(guī)有限元法的計(jì)算結(jié)果Table 2 Result of conventional finite element method

從表1、2可以看出，本文ODE解只用了3次迭代就求得彈性壓桿穩(wěn)定問題的精確解，而利用常規(guī)有限元法需要將結(jié)構(gòu)劃分為48個(gè)單元格，才求得相應(yīng)精度的解，說明對于相同精度的解答，本文方法的計(jì)算效率明顯高于常規(guī)有限元法。此外，該方法還具有使用方便、計(jì)算量少、收斂速度快等優(yōu)點(diǎn)。

4 結(jié)論

本文提出的常微分方程求解器解法，采用數(shù)次迭代就可以得到所求問題的精確解，不需要對單元進(jìn)行細(xì)分，而且該方法從計(jì)算精度和計(jì)算效率上均優(yōu)于常規(guī)有限元法，并且其精度可由用戶所任意指定，使用操作也相當(dāng)方便。由此可預(yù)見，這種常微分方程求解器數(shù)值解法的應(yīng)用前景是相當(dāng)廣闊的。

［1］楚中毅，陸念力，楚蘭英，等.梁桿結(jié)構(gòu)穩(wěn)定性分析的一種精確有限元方法及其優(yōu)化［J］.建筑機(jī)械，2001，9(14):68－72

［2］楚中毅，陸念力，車仁煒，等.一種梁桿結(jié)構(gòu)穩(wěn)定性分析的精確有限元法［J］.哈爾濱建筑大學(xué)學(xué)報(bào)，2002，4(5):25－28

［3］龍馭球，包世華.結(jié)構(gòu)力學(xué)［M］.北京:高等教育出版社，1999

［4］任鳳鳴，范學(xué)明.彈性壓桿穩(wěn)定問題的精確解法［J］.建筑科學(xué)，2008，24(3):12－14

［5］包世華.結(jié)構(gòu)力學(xué)Ⅱ［M］.北京:高等教育出版社，2001

［6］包世華，周 堅(jiān).薄壁桿件結(jié)構(gòu)力學(xué)［M］.北京:中國建筑出版社，1991

［7］A scher U，Christiansen J and Russell R D.Collocation Software for Boundary－value ODE［J］.ACM Trans Math Software，1981，7(2):209－222

［8］A scher U，Christiansen J and Russell R D.Algorithm 569，COLSYS:Collocation Sof tware for Boundary－value ODEs［J］.ACM Trans Math Sof tware，1981，7(2):223－229