基于FORTRAN的3D等效結點荷載計算

司馬丹琪++李元松++姜成潼++何泉

摘要： 在有限單元分析中，需將作用于單元的外力按虛功等效原則，移置到結點上成為等效結點力。由于數學運算的困難，往往限制了工程技術人員運用有限單元技術分析實際具有分布力的工程問題。以常用空間20結點單元為例，詳述荷載等效過程，附注自編FORTRAN程序。并利用FEA軟件進行計算對比驗證程序正確性，為學者學習溫度、預應力等類似工程問題的等效結點荷載計算提供有益參考。

Abstract： Abstract： In the finite element analysis， the external force acting on the element needs to be transferred to the node as the equivalent node force according to the principle of virtual work equivalence. Due to the difficulty of mathematical operations， engineers and technicians are often limited to use finite element technique to analyze practical distributed engineering problems. Taking the commonly used space 20-node element as an example， the equivalent process of load is described in detail. The software FEA is used to calculate and compare the correctness of the program， which provides a useful reference for scholars to study the equivalent nodal load of temperature， prestress and other similar engineering problems.

關鍵詞： 離散體；等效結點荷載；虛功原理；程序設計

Key words： discrete body；equivalent node load；virtual work principle；program design

中圖分類號：TP311.1 文獻標識碼：A 文章編號：1006-4311（2018）08-0233-02

0 引言

在彈性介質靜力問題的計算中外加作用因素很多，如集中力、分布力（包括引力、斥力等場力、慣性力一類體積力和面積分布的由力邊界條件給定的接觸力）等直接載荷，也可能有因為溫度改變、裝配因素、預應力作用等其它干擾力[1]。這些外加作用因素都可遵循力學等效原則（如靜力等效、位移模式下的虛功等效）處理成結點載荷。這種等效處理往往涉及較為復雜的坐標變換運算[2]，對于非力學專業(yè)的工程技術人員而言，存在一定困難。并且由于計算機語言的限制，傳統(tǒng)的有限元程序設計課程只能以二維有限元問題為例介紹程序設計過程，極少涉及三維有限元的編程。而實際工程問題無一不是三維問題，因此編寫三維程序更具有實際價值?，F有商業(yè)計算軟件中邊界處理條件功能均很強大，但很難涵蓋工程實際中遇到的各種邊界條件問題，一旦遇到軟件中沒有對應的處理方法，仍需根據工程條件自行開發(fā)程序解決。筆者在有限元程序系統(tǒng)（FEM-PS）開發(fā)過程中，對3D等參單元等效結點荷載的計算公式進行詳細推導，并利用FORTRAN程序設計成通用模塊，以期對類似問題的解決有所啟示。本文以三維空間二十結點等參單元為例，介紹FORTRAN語言在程序設計中的使用，并針對三維空間設計過程中的等效結點荷載計算環(huán)節(jié)做程序詳述。

1 等效結點荷載計算原理

1.1 集中力

設在單元內部或邊界上任意點C（x，y，z）作用有集中荷載{Q}=[Qx，Qy，Qz]T，轉化成等效結點荷載列陣{PE}（e）。根據轉換前后虛功相等原則有[3]：

{PE} （e）=[Nc]T{Q}（1）

式中[Nc]為形函數在C點的值。

1.2 體積力

設單元內作用有均布體積力{p}=[px，py，pz]T，則可將單元微分體積dV上的體力{p}dV視為集中荷載，利用式（1），在單元的體積V內積分，即得到該分布體積力的等效結點荷載{PE }

1.3 表面力

設在單元的某邊界面上作用表面力{ }={ ， ， }T，可將微分面積dA上的面荷載{p}dA視為集中荷載，利用式（1）在面積A上積分，即得到該分布面荷載的等效結點荷載{PE}

式中Ni表示單元的形函數，對于空間六面體20結點單元，其表達式為：

由于工程實際中給出的表面力往往是沿作用面的法向力和切向力，對于空間問題，外力往往垂直作用于邊界表面，以？孜=1的面為例，載荷集度以 0表示，{PE} （e）表達式變?yōu)閇4]

在計算結點荷載向量時由于被積函數比較復雜，常采用三維高斯積分公式[5]求解：

2 FORTRAN的程序實現

作用在單元上的面力在等效成結點荷載的過程中，由于坐標軸的確定，力所作用面的方向會發(fā)生相應的變化。本文對三維二十結點的六個空間面進行編號，新定義變量NOACE（*）來表示力作用的單元面；將力作用面固定，如EXISP=-1.0（ξ=-1.0），再利用FORTRAN程序內部存儲各面形函數順序數組，如NFACE=（/1，9，5，20，8，12，4，16/），確定空間六面體結點拓撲關系：

IF（NOACE.EQ.1） THEN ！面力作用單元面的選擇

NFACE=（/1，9，5，20，8，12，4，16/）

EXISP= -1.0

最后調用其他公共子模塊矩陣[6]，求解分布面力的等效結點荷載。利用斷點法對程序進行調試過程中，單元在受均布面力作用下進行等效之后，四周角點的等效結點荷載變?yōu)?1/12q，邊中點變?yōu)?/3q，所述規(guī)律與文獻[7]中相符。面力等效程序如下：

LGASH=（NFACE（IODEG）-1）*NDOFN+1 ！確定結點自由度編碼

MGASH=（NFACE（IODEG）-1）*NDOFN+2

NGASH=（NFACE（IODEG）-1）*NDOFN+3 ！荷載等效

ELOAD（NEASS，LGASH）=ELOAD（NEASS，LGASH）+SHAPEI（NFACE（IODEG））*PXCOM* WEIGP（IGAUS）*WEIGP（JGAUS）

ELOAD（NEASS，MGASH）=ELOAD（NEASS，MGASH）+SHAPEI（NFACE（IODEG））*PYCOM*WEIGP（IGAUS）*WEIGP（JGAUS）

ELOAD（NEASS，NGASH）=ELOAD（NEASS，NGASH）+SHAPEI（NFACE（IODEG））*PZCOM* WEIGP（IGAUS）*WEIGP（JGAUS）

3 實例及分析

某空間正方體壓力盒各邊長為2.0m，上端面受均布面力荷載q=10.0kN/m2，其余面均固定，E=2.0×105MPa，μ=0.2，體密度ρ=50.0kN/m3，采用20節(jié)點等參單元，單元個數8，節(jié)點總數81。（圖1）

為驗證程序的正確性，在劃分相同網格的條件下，取有限元軟件FEA計算結果與程序計算結果做比對，如圖2、圖3所示。

由圖2、圖3的對比分析可知，用FORTRAN編制的程序計算結果與利用FEA模擬出來的結果基本吻合，說明了編制程序的正確性。

4 結語

有限元商業(yè)軟件的核心仍舊是有限元算法，只有掌握相應的邏輯流程及算法原理才能熟練操作一款軟件，解決遇到的棘手問題。文章基于虛功原理推導了彈性問題空間20結點單元等效結點荷載計算公式，利用FORTRAN語言，設計成通用模塊程序。采用算例驗證算法與程序的正確性，為類似連續(xù)介質體問題的有限元法求解提供有益的參考。

參考文獻：

[1]徐芝綸.彈性力學簡明教程[M].三版.北京：高等教育出版社，1983.

[2]王勖成.有限單元法[M].北京：清華大學出版社，2003.

[3]趙更新，等.土木工程結構分析程序設計[M].北京：中國水利水電出版社，2001.

[4]王元漢，李麗娟，李銀平，等.有限元基礎與程序設計[M].廣州：華南理工大學出版社，2002.

[5]朱加銘，歐貴寶，何蘊增.有限元與邊界元法[M].哈爾濱：哈爾濱工程大學出版社，2002.

[6]Hinton E and D R J. Finite Element Programming[M].London： Academic Press，1977.

[7]朱伯芳.有限單元法原理與應用[M].二版.北京：中國水利水電出版社，1998.