非線性方程組求解器及平面連桿機(jī)構(gòu)仿真程序研究*

柳冬玉，潘亞嘉

(西南交通大學(xué) 機(jī)械工程學(xué)院，四川 成都 610031)

0 引 言

傳統(tǒng)的機(jī)構(gòu)設(shè)計(jì)有幾何法和解析法。在對(duì)某一特定機(jī)構(gòu)進(jìn)行參數(shù)分析時(shí)，可使用CAD軟件建模，進(jìn)而進(jìn)行機(jī)構(gòu)分析。為分析機(jī)構(gòu)而進(jìn)行三維建模耗工耗時(shí)，部分CAD軟件雖然具備骨架功能，但操作仍較繁瑣。對(duì)于機(jī)構(gòu)的運(yùn)動(dòng)學(xué)問題，手工計(jì)算或者編程求解又面臨求解非線性方程組的問題，研究簡明易用的平面連桿機(jī)構(gòu)仿真程序十分必要。

目前平面機(jī)構(gòu)仿真程序有以下幾類：

(1) 專門針對(duì)特定機(jī)構(gòu)的程序，首先讓用戶從幾種固定的機(jī)構(gòu)類型中做出選擇，再輸入尺寸進(jìn)行仿真，或者輸入設(shè)計(jì)目標(biāo)得到機(jī)構(gòu)尺寸。此類程序算法固定，不能夠處理內(nèi)置機(jī)構(gòu)以外的類型。

(2) 以桿組理論為核心的通用平面機(jī)構(gòu)仿真程序，用戶必須以桿組為基本單元進(jìn)行輸入[1-5]。

(3) 不依賴桿組，以單個(gè)桿件為基本元素的通用平面機(jī)構(gòu)仿真程序。用戶需輸入桿件及其約束信息[6-7]。

基于桿組理論的機(jī)構(gòu)仿真程序要求用戶熟悉桿組分解，并且不能夠求解程序中沒有的桿組，具有先天缺陷。筆者在上述第3類程序的基礎(chǔ)上，引入符號(hào)計(jì)算思想所進(jìn)行的研究。

1 基于約束的運(yùn)動(dòng)學(xué)模型建立

1.1 機(jī)構(gòu)的約束方程組的一般形式

設(shè)機(jī)構(gòu)由一個(gè)約束方程組描述[8]為：

(1)

由式(1)可以解出機(jī)構(gòu)在任意時(shí)間下的q，從而得出機(jī)構(gòu)的位置。由于式(1)通常為高度非線性，所以大多使用數(shù)值方法求解。

由于解出的q非解析解，所以不能直接對(duì)q微分得到機(jī)構(gòu)各構(gòu)件的速度及加速度。

(2)

(3)

考慮由轉(zhuǎn)動(dòng)副和移動(dòng)副構(gòu)成的平面連桿機(jī)構(gòu)，其基本約束形式為重合約束和共線約束。重合約束的方程形式為：

(4)

式中：xi,yi,φi為構(gòu)件i的笛卡爾廣義坐標(biāo)；xj,yj,φj為構(gòu)件j的笛卡爾廣義坐標(biāo)；xi′P,yi′P,xj′P,yj′P為兩構(gòu)件重合點(diǎn)P分別在構(gòu)件i,j的固定坐標(biāo)系上的坐標(biāo)。若取其中一個(gè)構(gòu)件的廣義坐標(biāo)為常數(shù)，則可得到機(jī)架點(diǎn)鉸接活動(dòng)構(gòu)件的約束方程。

共線約束的方程形式為：

xi′Q)]×(xi-xj+xi′Pcosφi-xj′Pcosφj-

yi′Psinφi+yj′Psinφj)-[cosφi(xi′P-

xi′Q)-sinφi(yi′P-yi′Q)]×(yi-yj+

yi′Pcosφi-yj′Pcosφj+xi′Psinφi-

xj′Psinφj)≡0

Φ2t(i,j)=[cos(φi-φj)(yi′P-yi′Q)+

sin(φi-φj)(xi′P-xi′Q)](xj′P-xj′Q)-

[cos(φi-φj)(xi′P-xi′Q)-sin(φi-

(5)

式中：xi,yi,φi,xj,yj,φj含義同式(4)；兩構(gòu)件間存在一條公共的線段PQ，點(diǎn)P和點(diǎn)Q在構(gòu)件i,j的固定坐標(biāo)系上的坐標(biāo)分別為xi′P,yi′P,xj′P,yj′P和xi′Q,yi′Q,xj′Q,yj′Q。若取其中一個(gè)構(gòu)件的廣義坐標(biāo)為常數(shù)，則可得到機(jī)架與活動(dòng)構(gòu)件共線的方程。

1.2 拖動(dòng)約束方程

很多CAD軟件都具有鼠標(biāo)拖動(dòng)零件查看機(jī)構(gòu)運(yùn)動(dòng)的功能，拖動(dòng)功能能夠非常直觀地反映機(jī)構(gòu)的運(yùn)動(dòng)情況，允許用戶拖動(dòng)任一構(gòu)件查看仿真效果。下面推導(dǎo)拖動(dòng)約束方程。如圖1所示，構(gòu)件i在點(diǎn)MStart處被鼠標(biāo)點(diǎn)擊拖動(dòng)至MDrag，構(gòu)件受到其余約束方程的限制進(jìn)行運(yùn)動(dòng)，廣義坐標(biāo)發(fā)生變化，但運(yùn)動(dòng)規(guī)律未知。

圖1 構(gòu)件i拖動(dòng)前后示意

發(fā)生拖動(dòng)前后，點(diǎn)M與構(gòu)件的相對(duì)位置不發(fā)生變化，點(diǎn)M在坐標(biāo)系xi′-yi′上的坐標(biāo)為:

SM′=[xm′ym′]T

(6)

矢量SM與坐標(biāo)軸xi′的夾角為:

(7)

點(diǎn)MDrag存在關(guān)系:

l=(xm-xi)sin(φi+α)

=(ym-yi)cos(φi+α)

(8)

式中：l為點(diǎn)M至坐標(biāo)系xi′-yi′距離。

整理，得到拖動(dòng)約束方程為:

(9)

已知一個(gè)機(jī)構(gòu)系統(tǒng)的所有尺寸及其初始狀態(tài)，列出所有轉(zhuǎn)動(dòng)副與移動(dòng)副的約束方程，與驅(qū)動(dòng)約束方程聯(lián)立可得式(1)，依次求解式(1)～(3)則可得到機(jī)構(gòu)系統(tǒng)任意時(shí)刻的位移、速度及加速度。若與拖動(dòng)約束方程聯(lián)立求解，則可得到機(jī)構(gòu)系統(tǒng)在受約束限制時(shí)拖動(dòng)某一構(gòu)件后的位置。

2 利用表達(dá)式樹構(gòu)建非線性約束方程組求解器

由于運(yùn)動(dòng)學(xué)方程是高度非線性的，其求解方式至關(guān)重要。基于牛頓—拉夫森迭代法的運(yùn)動(dòng)學(xué)約束方程求解在程序編制上通常使用固定矩陣法，其步驟如下：

(1) 事先展開所有運(yùn)動(dòng)副及驅(qū)動(dòng)約束的位移、速度、加速度方程，推導(dǎo)特定約束對(duì)于兩個(gè)構(gòu)件的位移、速度、加速度方程右項(xiàng)。

(2) 在計(jì)算時(shí)代入構(gòu)件的廣義坐標(biāo)，組合成位移、速度、加速度整體的數(shù)值矩陣，求解線性方程組后再根據(jù)迭代終止條件回代到位移、速度、加速度右項(xiàng)中繼續(xù)迭代。

固定矩陣法整個(gè)求解過程不涉及符號(hào)推導(dǎo)，運(yùn)算速度快，但是因?yàn)槊糠N方程的計(jì)算表達(dá)式已經(jīng)推導(dǎo)并且固化在程序中，所以無法為原動(dòng)件設(shè)置自定義函數(shù)的運(yùn)動(dòng)規(guī)律，也無法求解運(yùn)動(dòng)學(xué)方程以外的方程。文獻(xiàn)[8]討論的即是此類求解器。

鑒于固定矩陣法的缺陷，引入符號(hào)計(jì)算思想及編譯原理編制了通用非線性方程組求解器核心，在求解器核心的基礎(chǔ)上添加自動(dòng)列式功能及接口構(gòu)成非線性約束方程組求解器。求解器可識(shí)別中序遍歷表達(dá)式的非線性方程組字符串，并求出結(jié)果。文中的方法與固定矩陣法的求解流程圖對(duì)比如圖2所示。

圖2 固定矩陣法與文中方法求解流程對(duì)比

求解器功能由四個(gè)類實(shí)現(xiàn)：求解器類(TSolver)、方程組類(TEquations)、表達(dá)式樹類(TExpressionTree)、變量表類(TVariableTable)。

變量表類用于保存變量名及其初始值；表達(dá)式樹類用于將方程字符串解析為表達(dá)式樹；方程組類用于存儲(chǔ)各個(gè)方程的表達(dá)式樹并對(duì)其求解，可獨(dú)立求解非線性方程組；求解器類用于讀取程序構(gòu)件數(shù)據(jù)建立約束方程，對(duì)外提供方程組類中的求解結(jié)果的接口。

在方程組類內(nèi)部，方程字符串由詞法分析器解析為表達(dá)式樹，表達(dá)式樹進(jìn)行符號(hào)求導(dǎo)得到由表達(dá)式樹組成的雅可比(Jacobian)符號(hào)矩陣。雅可比符號(hào)矩陣及方程組均代入初值得到線性方程組，利用高斯消元中的列主元消元法得到Δq，同時(shí)判斷出雅可比矩陣的奇異性。圖3為基于表達(dá)式樹的非線性方程組求解過程圖。

圖3 基于表達(dá)式樹的非線性方程組求解過程

根據(jù)Δq進(jìn)行收斂判斷，若Δq大于設(shè)定值則繼續(xù)迭代。若Δq小于設(shè)定精度則終止迭代并輸出求解結(jié)果。

3 程序結(jié)構(gòu)組成及設(shè)計(jì)

3.1 程序結(jié)構(gòu)

程序以面向?qū)ο竽Ｊ介_發(fā)，其結(jié)構(gòu)分為界面層、工具層、構(gòu)件層和公共層，如圖4所示。

圖4 程序結(jié)構(gòu)組成

界面層負(fù)責(zé)程序的界面顯示，各個(gè)控件的業(yè)務(wù)邏輯，彈出對(duì)話框。工具層包含工具庫，繪制類，吸附類和求解器：工具庫包含實(shí)現(xiàn)各類構(gòu)件繪制的工具類；繪制類提供各類構(gòu)件及矢量圖形的繪制函數(shù)；吸附類用于在構(gòu)件繪制時(shí)提供吸附功能；求解器如前文所述，用于求解運(yùn)動(dòng)學(xué)約束方程。構(gòu)件層將各類構(gòu)件、約束及驅(qū)動(dòng)封裝為類，構(gòu)件管理器用于管理構(gòu)件信息，輸入輸出構(gòu)件數(shù)據(jù)文件。公共層包含配置類、公共函數(shù)庫及資源，配置類存儲(chǔ)程序當(dāng)前的單位、比例、坐標(biāo)原點(diǎn)、線型、顏色等信息，供其他層調(diào)用。

界面層與工具層在實(shí)現(xiàn)邏輯上通過直接調(diào)用Windows API實(shí)現(xiàn)，不依賴JRE、.net Framework等框架及MFC、Qt等類庫。

3.2 求解算法及動(dòng)畫邏輯

構(gòu)建完成機(jī)構(gòu)系統(tǒng)后，拖動(dòng)構(gòu)件、點(diǎn)擊刷新、運(yùn)行仿真都將激活求解器。求解時(shí)，程序?qū)?gòu)件、約束及驅(qū)動(dòng)進(jìn)行遍歷，生成約束方程及驅(qū)動(dòng)方程。此時(shí)得到整個(gè)機(jī)構(gòu)的方程組，之后將機(jī)構(gòu)的當(dāng)前位置作為求解器中各未知量的初值，代入求解器進(jìn)行求解。若求解成功，則將求解器中各未知量數(shù)據(jù)傳回桿件，刷新畫布顯示結(jié)果。程序支持動(dòng)畫功能。將時(shí)間參數(shù)連續(xù)傳入求解器并存儲(chǔ)求解結(jié)果，將求解結(jié)果連續(xù)送入畫布繪制機(jī)構(gòu)，即可實(shí)現(xiàn)動(dòng)畫效果。

4 應(yīng)用實(shí)例

牛頭刨床六桿機(jī)構(gòu)如圖5所示。

圖5 牛頭刨床六桿機(jī)構(gòu)運(yùn)動(dòng)簡圖

各構(gòu)件尺寸為l1=125 mm，l3=600 mm，l4=150 mm，主動(dòng)件為構(gòu)件1，以ω1=360 (°)/s等速轉(zhuǎn)動(dòng)，初始位置為φ1=0 (°)，l5為從動(dòng)件，求E點(diǎn)處x方向的位移、速度及加速度。

建模及分析過程如下。首先在圖6所示程序主界面上點(diǎn)選工具，在操作區(qū)繪制各構(gòu)件，繪制過程中將自動(dòng)產(chǎn)生捕捉，并建立起約束關(guān)系，操作方式類似于AutoCAD。

圖6 程序界面

之后使用選擇工具選中連桿1，點(diǎn)擊工具欄上“設(shè)為原動(dòng)件”按鈕，彈出如圖7所示窗口，設(shè)置原動(dòng)件運(yùn)動(dòng)規(guī)律。此時(shí)機(jī)構(gòu)建立完畢。

圖7 設(shè)置原動(dòng)件界面

點(diǎn)擊主界面工具欄上“動(dòng)畫與測量”按鈕，彈出窗口如圖8所示。設(shè)置好各測量項(xiàng)目，點(diǎn)擊“開始分析”按鈕，程序?qū)?duì)機(jī)構(gòu)在指定時(shí)間內(nèi)的運(yùn)動(dòng)進(jìn)行分析。分析完成后，將彈出測量結(jié)果線圖，即E點(diǎn)在x方向的位移、速度及加速度，如圖9所示。點(diǎn)擊線圖窗口菜單項(xiàng)可將線圖保存為圖片及csv數(shù)據(jù)文件。在“動(dòng)畫與測量”界面上拖動(dòng)時(shí)間軸或者點(diǎn)擊播放等按鈕可以查看機(jī)構(gòu)運(yùn)動(dòng)仿真動(dòng)畫，如圖10。

圖8 動(dòng)畫與測量界面

圖9 LML輸出位移、速度、加速度線圖

圖10 ADAMS輸出位移、速度、加速度線圖

表1為本文程序LML進(jìn)行位移、速度、加速度分析的數(shù)據(jù)與ADMAS的對(duì)比。從線圖及輸出數(shù)據(jù)來看，牛頭刨床推程運(yùn)動(dòng)速度平穩(wěn)，加速度小；回程速度變化大，加速度大，符合急回運(yùn)動(dòng)特性，驗(yàn)證了程序求解的正確性。

表1 位移、速度、加速度分析數(shù)據(jù)與ADAMS對(duì)比

5 結(jié) 論

建立了基于約束的平面連桿機(jī)構(gòu)系統(tǒng)運(yùn)動(dòng)學(xué)模型，通過表達(dá)式樹實(shí)現(xiàn)基于牛頓—拉夫森迭代法的通用非線性方程組求解器，開發(fā)了平面連桿機(jī)構(gòu)仿真程序。

(1) 封裝了Windows基本窗口及原生控件，內(nèi)部功能完全通過Windows API實(shí)現(xiàn)，程序平臺(tái)依賴性低，體積小，可在沒有安裝任何運(yùn)行庫的Windows XP以上系統(tǒng)運(yùn)行。封裝的窗口類及控件類可作為獨(dú)立類庫用于開發(fā)Win32界面程序。

(2) 求解器核心基于C++面向?qū)ο缶幊谭椒ㄔO(shè)計(jì)，可作為通用的線性、非線性及不定方程組解算工具。求解器不依賴桿組理論，可求解任意桿組級(jí)別、多自由度機(jī)構(gòu)的運(yùn)動(dòng)學(xué)約束方程。

(3) 圍繞平面連桿機(jī)構(gòu)的設(shè)計(jì)、編輯、分析、動(dòng)畫功能設(shè)計(jì)了簡明的用戶操作界面及人機(jī)交互方式，實(shí)現(xiàn)真正的“拖拽式編輯”，程序自動(dòng)吸附節(jié)點(diǎn)、添加約束，用戶不需要介入節(jié)點(diǎn)編號(hào)等繁瑣細(xì)節(jié)，所有操作所見即所得，并可拖動(dòng)機(jī)構(gòu)構(gòu)件查看任意運(yùn)行位置，原動(dòng)件可設(shè)置任意表達(dá)式的運(yùn)動(dòng)規(guī)律，輸出位移、速度、加速度圖像及數(shù)據(jù)。與ADAMS的仿真結(jié)果對(duì)比驗(yàn)證了求解的正確性。

[1] 李永旭,劉 釗,李小江.平面連桿機(jī)構(gòu)計(jì)算機(jī)建模與通用分析軟件的開發(fā)[J].機(jī)械設(shè)計(jì),2006,23(5): 58-60.

[2] 田 松.基于桿組理論的III級(jí)機(jī)構(gòu)運(yùn)動(dòng)分析與仿真[D].武漢:武漢科技大學(xué),2007.

[3] 李小江.MAD中的RRR桿組原理及應(yīng)用[J].機(jī)械設(shè)計(jì)與制造,2008(3):66-68.

[4] 黃 河.平面連桿機(jī)構(gòu)通用分析程序的研究與開發(fā)[D].西安: 長安大學(xué),2008.

[5] 岳一領(lǐng),陸鳳儀,張志鴻.連桿機(jī)構(gòu)設(shè)計(jì)仿真平臺(tái)在機(jī)械原理教學(xué)中的應(yīng)用[J].太原理工大學(xué)學(xué)報(bào)(社會(huì)科學(xué)版),2009,27(1):74-76, 88.

[6] 王炎歡.機(jī)構(gòu)仿真軟件的研究與開發(fā)[D].合肥:合肥工業(yè)大學(xué),2006.

[7] 徐春雷，潘宏俠，李 冉，等.基于Matlab的平面連桿機(jī)構(gòu)仿真軟件的研究[J].機(jī)械設(shè)計(jì)與制造,2012,12(12):94-96.

[8] Haug Edward J. Computer aided kinematics and dynamics of mechanical systems:Basic methods Vol.1[M].Boston:Allyn and Bacon, 1989.

[9] Charles Petzold. Programming Windows: 5thed.[M]. Redmond: Microsoft Press,1998.

[10] Stanley B. Lippman, Josée Lajoie, Barbara E. Moo. C++ Primer: 5thed[M].Boston: Addison-Wesley Professional, 2012.

[11] Weifeng Huang. Automated synthesis of planar mechanisms with revolute, prismatic and pin-in-slot joints[D]. United States: Oregon State University, 2015.

[12] Feng Yuan. Windows graphics programming: Win32 GDI and DirectDraw[M].New Jersey: Prentice Hall, 2000.

[13] 李正彪,余 波.用矩陣的初等列變換求解多元線性不定方程[J].云南民族大學(xué)學(xué)報(bào)(自然科學(xué)版),2006,15(2):102-104.

[14] 李 超,阮 威,張 龍,等.計(jì)算機(jī)代數(shù)系統(tǒng)數(shù)學(xué)原理[M].北京:清華大學(xué)出版社,2010.

[15] 金一慶, 陳 越, 王冬梅.數(shù)值方法[M].第2版.北京:機(jī)械工業(yè)出版社,2015.