張明哲，彭妙娟，程玉民

(上海大學(xué) 力學(xué)與工程科學(xué)學(xué)院，上海 200444)

0 引 言

在求解高速碰撞、動(dòng)態(tài)裂紋擴(kuò)展和材料幾何非線性等問題時(shí)，無(wú)網(wǎng)格方法可以基于一系列離散的節(jié)點(diǎn)構(gòu)造試函數(shù)，從根本上避免采用有限元法容易發(fā)生的網(wǎng)格畸變、計(jì)算過(guò)程終止等缺點(diǎn)。無(wú)網(wǎng)格方法受到越來(lái)越多國(guó)內(nèi)外學(xué)者的關(guān)注，因此開發(fā)一款以開放式語(yǔ)言為基礎(chǔ)的無(wú)網(wǎng)格方法計(jì)算軟件很有必要。

有關(guān)彈塑性問題的無(wú)網(wǎng)格方法軟件開發(fā)起步較晚，現(xiàn)階段比較完善且運(yùn)用廣泛的工程分析軟件仍然是有限元法軟件，如MSC Nastran、Abaqus、Ansys、MIDAS等。Nayroles等基于移動(dòng)最小二乘法提出擴(kuò)散單元法，無(wú)網(wǎng)格方法才真正發(fā)展起來(lái)。王宇新等開發(fā)無(wú)網(wǎng)格MPM法三維前處理系統(tǒng)；張征等開發(fā)接觸力學(xué)自適應(yīng)無(wú)網(wǎng)格計(jì)算系統(tǒng)，求解二維彈塑性接觸力學(xué)；ABBASZADEH等開發(fā)一種無(wú)網(wǎng)格數(shù)值程序模擬地下水污染方程；曾清紅利用并行計(jì)算技術(shù)進(jìn)行適于無(wú)網(wǎng)格方法的并行插值算法研究；NAKATA等提出一種在圖形處理單元上的基于徑向點(diǎn)插值方法的無(wú)網(wǎng)格分析并行算法；SINGH等基于Fortran語(yǔ)言開發(fā)并行算法，實(shí)現(xiàn)傳熱和流體流動(dòng)問題的無(wú)單元Galerkin方法的并行求解；LIU等使用MLS為無(wú)網(wǎng)格方法開發(fā)一種基于背景網(wǎng)格的自適應(yīng)程序；文建波等開發(fā)基于Delaunay三角化的無(wú)網(wǎng)格方法計(jì)算結(jié)果后處理軟件。此外，Ansys作為成熟的有限元軟件，已經(jīng)內(nèi)置無(wú)網(wǎng)格計(jì)算方法的模塊。

關(guān)于彈塑性問題的無(wú)網(wǎng)格方法的研究，傳統(tǒng)無(wú)網(wǎng)格方法采用移動(dòng)最小二乘法構(gòu)造試函數(shù)，在計(jì)算過(guò)程中容易形成病態(tài)或者奇異方程組，而彈塑性問題的求解要求布置足夠多數(shù)量的節(jié)點(diǎn)以提高計(jì)算精度，這嚴(yán)重影響計(jì)算效率和精度。為解決這些問題并擴(kuò)大無(wú)網(wǎng)格方法的應(yīng)用范圍，CHEN等應(yīng)用復(fù)變量重構(gòu)核粒子法，用一維基函數(shù)形成二維問題的校正函數(shù)解決二維彈塑性問題；CHENG等使用可直接施加邊界條件的插值移動(dòng)最小二乘法獲得形函數(shù)，提出插值無(wú)單元Galerkin方法解決二維彈塑性問題；CHENG等采用罰函數(shù)法應(yīng)用本質(zhì)邊界條件解決二維彈塑性問題，獲得更高的計(jì)算精度和效率；MENG等利用維數(shù)分裂法將三維波動(dòng)方程分解為二維問題，推導(dǎo)改進(jìn)的插值維數(shù)分裂無(wú)單元Galerkin方法求解三維波動(dòng)方程；蔡小杰等采用移動(dòng)最小二乘法建立彈塑性大變形問題的改進(jìn)的無(wú)單元Galerkin方法；DENG等基于改進(jìn)的復(fù)變量無(wú)單元Galerkin方法提出二維非線性的彈塑性問題的數(shù)值模型；WU等、PENG等和ZHENG等提出改進(jìn)的無(wú)單元Galerkin方法解決彈塑性大變形問題，利用Galerkin積分弱形式建立求解方程，其形函數(shù)滿足Kronecker函數(shù)的性質(zhì)，可直接施加本質(zhì)邊界條件，并提出采用插值無(wú)單元Galerkin方法求解三維彈塑性問題，可極大地提高計(jì)算效率和精度，由其編寫的無(wú)網(wǎng)格方法部分程序算法，可極大地減少計(jì)算量。

上述關(guān)于彈塑性問題的無(wú)網(wǎng)格方法的程序大多不具有全流程分析(精確的計(jì)算運(yùn)行程序和可視化的前、后處理程序)，如MATLAB專注于計(jì)算程序，使得分析問題時(shí)不能在一個(gè)軟件中兼顧建模、計(jì)算分析和結(jié)果展示，也間接阻礙彈塑性問題無(wú)網(wǎng)格方法的普及和應(yīng)用。無(wú)網(wǎng)格方法的相關(guān)軟件很少有關(guān)于彈塑性分析的。

本文開發(fā)彈塑性問題的改進(jìn)無(wú)單元Galerkin方法軟件，包括前處理程序模塊、計(jì)算程序模塊和后處理程序模塊3個(gè)部分。前處理程序鏈接Gmsh建模，可以獲得節(jié)點(diǎn)信息等信息流，然后通過(guò)彈塑性問題的改進(jìn)無(wú)單元Galerkin方法計(jì)算程序進(jìn)行數(shù)值模擬，最后利用后處理程序模塊實(shí)現(xiàn)數(shù)值結(jié)果的可視化。

1 軟件架構(gòu)

該軟件基于C++編程語(yǔ)言開發(fā)，依托改進(jìn)的無(wú)單元Galerkin方法解決彈塑性問題，能夠高效直觀地建模并進(jìn)行準(zhǔn)確的受力分析。軟件面向的對(duì)象包括節(jié)點(diǎn)、單元以及邊界條件等模型數(shù)據(jù)，對(duì)其進(jìn)行計(jì)算分析后得到處理結(jié)果，最后將結(jié)果在軟件上以可視化界面進(jìn)行顯示，即前處理程序、分析計(jì)算程序和后處理程序3個(gè)部分。前處理程序、計(jì)算程序和后處理程序不是互相獨(dú)立工作的，其數(shù)據(jù)聯(lián)系流程見圖1。

圖1 軟件流程圖

前處理程序模塊主要為建立模型和節(jié)點(diǎn)生成：基于OpenGL庫(kù)直接使用C++語(yǔ)言進(jìn)行建模，輸入?yún)?shù)，預(yù)處理必備的結(jié)構(gòu)參數(shù)值。因無(wú)網(wǎng)格方法需要節(jié)點(diǎn)進(jìn)行函數(shù)擬合，其中的關(guān)鍵環(huán)節(jié)是對(duì)計(jì)算域離散剖分并獲取節(jié)點(diǎn)信息。無(wú)網(wǎng)格方法將計(jì)算域離散為一系列點(diǎn)，根據(jù)邊界法向前推進(jìn)，在計(jì)算域內(nèi)部進(jìn)行填充式布點(diǎn)，同時(shí)利用Delaunay算法進(jìn)行網(wǎng)格篩選以獲得節(jié)點(diǎn)信息。無(wú)網(wǎng)格方法分析計(jì)算部分首先根據(jù)建模的維度選擇二維彈塑性或三維彈塑性問題的改進(jìn)無(wú)單元Galerkin方法的程序，將計(jì)算得到的數(shù)據(jù)傳遞給后處理模塊，后處理模塊主要顯示結(jié)果。

前處理程序模塊需要輸入模型的幾何信息和物理信息，包括各種屬性和邊界參數(shù)等，建立分析模型后對(duì)模型進(jìn)行節(jié)點(diǎn)處理。前處理的重點(diǎn)是得到節(jié)點(diǎn)信息并讀入模型節(jié)點(diǎn)數(shù)據(jù)，為無(wú)網(wǎng)格方法分析做準(zhǔn)備。

本軟件鏈接簡(jiǎn)單的CAD引擎，通過(guò)編程先設(shè)定模型所需的主要節(jié)點(diǎn)，連點(diǎn)成線，然后生成面，最后生成體。每個(gè)點(diǎn)、線、面和體都有唯一的標(biāo)簽。程序可直接輸入各種幾何信息和物理信息，提供所需的各種屬性和邊界參數(shù)，如載荷大小以及坐標(biāo)、邊界條件等。

計(jì)算程序模塊進(jìn)行無(wú)網(wǎng)格方法的計(jì)算，最主要的是矩陣計(jì)算和求解。

后處理程序模塊對(duì)每個(gè)單元建立單元節(jié)點(diǎn)位移向量。以文件形式讀入前處理獲得的節(jié)點(diǎn)幾何信息和網(wǎng)格信息等，將計(jì)算程序得到的位移輸入到后處理的坐標(biāo)修正程序，通過(guò)改變節(jié)點(diǎn)和網(wǎng)格的坐標(biāo)值顯示模型的變形量，從而實(shí)現(xiàn)模型變形的可視化。

2 前后處理程序模塊

2.1 前處理程序

前處理是數(shù)值計(jì)算中非常重要的環(huán)節(jié)，不僅決定著程序解決問題的通用性，而且直接影響后續(xù)數(shù)值計(jì)算的精度和效率。前處理過(guò)程首要的是布置節(jié)點(diǎn)，包括節(jié)點(diǎn)編號(hào)及相應(yīng)的坐標(biāo)值。節(jié)點(diǎn)數(shù)據(jù)的空間位置信息定義為類P_Point，其C++程序如下。

class P_Point

{

public:

P_Point(const double x = 0,const double y = 0,

const double z = 0):_x(x),_y(y),_z(z) { }

～P_Point() { }

const double get X() const { return _x;}

const double get Y() const { return _y;}

const double get Z() const { return _z;}

void set X(const double x) { _x = x;}

void set Y(const double y) { _y = y;}

void set Z(const double z) { _z = z;}

bool operator <(const P_Point&p) const;

protected:

double *Point_p;

double _x,_y,_z;};

節(jié)點(diǎn)數(shù)據(jù)的獲取采用coord的vector容器接收。節(jié)點(diǎn)的初始三維空間信息存儲(chǔ)在vector容器中，每3個(gè)連續(xù)數(shù)據(jù)為1組。為方便使用，按照節(jié)點(diǎn)排布順序重新存儲(chǔ)數(shù)據(jù)，其coord的C++程序如下。

std::vectorcoord;

P_Point p;

for (int i = 0;i

{

double x = coord[i];

double y = coord[i + 1];

double z = coord[i + 2];

p = P_Point(x,y,z);

nodes.push_back(p);

}

前處理分為3個(gè)部分：使用C++編程建模；輸入材料的各要素?cái)?shù)值，如截面形狀及大小、載荷定義等；生成無(wú)網(wǎng)格方法分析所需要的節(jié)點(diǎn)數(shù)據(jù)。

在無(wú)網(wǎng)格方法前處理程序中，鏈接OpenGL庫(kù)直接使用C++語(yǔ)言建模，輸入各種幾何信息和物理信息，對(duì)模型信息屬性和邊界進(jìn)行配置，以便建模成功后生成網(wǎng)格。網(wǎng)格應(yīng)進(jìn)行合適性判斷，減少節(jié)點(diǎn)距離過(guò)近、避免出現(xiàn)丟失邊界及穿透邊界的情況，再通過(guò)邊交換方法對(duì)三角形進(jìn)行優(yōu)化，提高三角形質(zhì)量。

在求解域中選取適當(dāng)分布的離散節(jié)點(diǎn)，根據(jù)一定的準(zhǔn)則進(jìn)行篩選以形成適合計(jì)算的結(jié)構(gòu)。以二維問題為例，需要先確定一個(gè)點(diǎn)為中心點(diǎn)，再根據(jù)距離準(zhǔn)則篩選出一定數(shù)量的衛(wèi)星點(diǎn)，在這些衛(wèi)星點(diǎn)中運(yùn)用Delaunay準(zhǔn)則進(jìn)行臨時(shí)三角化處理，遍歷找出最合適的點(diǎn)云結(jié)構(gòu)，最后采用改進(jìn)的移動(dòng)最小二乘法進(jìn)行節(jié)點(diǎn)數(shù)據(jù)處理。

2.2 后處理程序

后處理程序讀取前處理的節(jié)點(diǎn)、單元等信息，獲取計(jì)算程序得到的位移等信息，并重新組織成標(biāo)準(zhǔn)后處理結(jié)果文件的形式，利用計(jì)算機(jī)可視化技術(shù)再現(xiàn)模型信息，將分析結(jié)果以圖像形式顯示出來(lái)。

本系統(tǒng)的后處理程序包括數(shù)據(jù)處理程序、圖像顯示程序等，具體內(nèi)容包括3個(gè)部分：(1)模型數(shù)據(jù)(初始的節(jié)點(diǎn)編號(hào)、坐標(biāo)等信息，以及計(jì)算程序得到的數(shù)據(jù))的提??；(2)模型結(jié)果可視化處理，分析結(jié)果數(shù)據(jù)，重新構(gòu)建模型的計(jì)算結(jié)果視圖；(3)結(jié)果數(shù)據(jù)查看，利用Visual Studio軟件的文件流函數(shù)將生成的數(shù)據(jù)直接存入文本文件。

前處理和計(jì)算程序產(chǎn)生的結(jié)果數(shù)據(jù)通常包括模型數(shù)據(jù)和結(jié)果數(shù)據(jù)，后處理程序采用文本文件存儲(chǔ)以上信息。文件接收網(wǎng)格節(jié)點(diǎn)等信息，并將其存儲(chǔ)到文本文件中，同時(shí)改變節(jié)點(diǎn)的空間坐標(biāo)以達(dá)到將變形可視化的目標(biāo)。文本文件網(wǎng)格節(jié)點(diǎn)等信息修正、顯示位移的操作程序如下。

//進(jìn)行文件操作

int _total = getFileLine("*****");//獲取文件

總行數(shù)函數(shù)

int Nodes_row=getStrrow("*****","*****");

//獲取某文件的某一指定字符串的行數(shù)

ifstream fin("*****");

if (!fin) // 如果讀取失敗，打印fail

{

cerr <<"fail" <

return;

}

ofstream fout("*****.txt");//創(chuàng)建一個(gè)新文件接收字符串

string s;

int _index = 0;

float data;

while (getline(fin,s)) {

if (_index <= Nodes_row || _index >=

EndNodes_row)

fout <

else{

istringstream iss(s);

vectornewdata;

while (iss >>data)

{

newdata.push_back(data);

}

bool is_chuan = false;

int _L = newdata.size();

if (_L == 3){

double chuandi = dataprocess(newdata[0],

deformedpoint);

fout <

else fout <

_index++;}

3 計(jì)算程序

3.1 改進(jìn)的移動(dòng)最小二乘法

采用基于點(diǎn)的近似構(gòu)造逼近函數(shù)是無(wú)網(wǎng)格方法的關(guān)鍵。形成逼近函數(shù)的主要方法有光滑粒子法、移動(dòng)最小二乘法和重構(gòu)核粒子法等，其中移動(dòng)最小二乘法應(yīng)用最為廣泛，即

(1)

式中：為基函數(shù)的個(gè)數(shù)；()為基函數(shù)；()為待定系數(shù)。線性基函數(shù)

=[1]

(2)

點(diǎn)鄰域內(nèi)的局部逼近為

(3)

定義泛函數(shù)

(4)

式中：為點(diǎn)影響域內(nèi)的節(jié)點(diǎn)；為影響域覆蓋點(diǎn)的節(jié)點(diǎn)數(shù)；(-)為定義的權(quán)函數(shù)。

()=()()

(5)

其中：

()=()

(6)

()=()

(7)

由以上計(jì)算可得逼近函數(shù)

(8)

式中：()為形函數(shù)，

()=[()() …()]=

()()()

(9)

在移動(dòng)最小二乘法中，當(dāng)較大時(shí)，式(5)有時(shí)是病態(tài)的甚至是奇異的，難以獲得正確解或者難以求解。由此，采用改進(jìn)的移動(dòng)最小二乘法，選取正交函數(shù)作為基函數(shù)，得到的方程不但可以避免病態(tài)和奇異，而且不需要求解方程的逆，可直接得到該方程組的解

()=()(),

(10)

其中：

(11)

將式(10)代入式(3)，可得

(12)

(13)

改進(jìn)的移動(dòng)最小二乘法求解系數(shù)()時(shí)不需要求解矩陣()的逆，可避免出現(xiàn)的病態(tài)或者奇異的方程組，因此可提高計(jì)算效率和精度。

3.2 彈塑性問題的改進(jìn)無(wú)單元Galerkin方法

基于改進(jìn)的移動(dòng)最小二乘法推導(dǎo)出二維彈性大變形問題和三維彈塑性問題的改進(jìn)無(wú)單元Galerkin方法的方程。三維問題可囊括二維問題相關(guān)公式，對(duì)三維公式推導(dǎo)的介紹更具有代表性，所以此處只介紹三維彈塑性問題的方程推導(dǎo)過(guò)程。

彈塑性問題的平衡方程可表示為

(Δ)+Δ=0,∈,

(14)

式中：(·)為微分算子矩陣；Δ為域內(nèi)任一點(diǎn)的應(yīng)力增量；Δ為域內(nèi)任一點(diǎn)的體力增量。

(15)

Δ=[ΔΔΔΔΔΔ]

(16)

Δ=[ΔΔΔ]

(17)

彈塑性問題的幾何方程表示為

Δ=(Δ)

(18)

式中：Δ為影響域內(nèi)任一點(diǎn)的應(yīng)變?cè)隽?；?span id="syggg00" class="emphasis_italic">為影響域內(nèi)任一點(diǎn)的位移增量。

彈塑性問題的物理方程表示為

Δ=Δ

(19)

在彈性階段內(nèi)和塑性階段內(nèi)表達(dá)式不同。

定義邊界條件為

(20)

(21)

在采用改進(jìn)的無(wú)單元Galerkin方法求解彈塑性問題時(shí)，需要使用適當(dāng)?shù)姆椒▽⑶蠼鈫栴}線性化，即使用一系列的線性解逼近非線性問題的解，例如采用逐步增加載荷的方法，按比例施加載荷。

由改進(jìn)的移動(dòng)最小二乘法的形函數(shù)表達(dá)式可得，域內(nèi)任一點(diǎn)的位移增量向量可表示為

(22)

式中：為點(diǎn)影響域內(nèi)的節(jié)點(diǎn)數(shù)。

應(yīng)變張量可表示為

(23)

綜合上述推導(dǎo)，可得

(24)

對(duì)表達(dá)式進(jìn)行替換，可得

(25)

(26)

(27)

(28)

(29)

若在應(yīng)力邊界上點(diǎn)處作用集中力，則

(30)

式(26)可轉(zhuǎn)化為

(31)

3.3 彈塑性問題的改進(jìn)無(wú)單元Galerkin方法計(jì)算程序開發(fā)

3.3.1 數(shù)據(jù)處理

模型進(jìn)行數(shù)值計(jì)算一般需要3類數(shù)據(jù)，分別為節(jié)點(diǎn)數(shù)據(jù)、單元數(shù)據(jù)和邊界條件數(shù)據(jù)。

節(jié)點(diǎn)數(shù)據(jù)主要包括節(jié)點(diǎn)的空間信息、編號(hào)，以及其對(duì)應(yīng)的計(jì)算得到的位移等數(shù)據(jù)。

單元數(shù)據(jù)主要包括定義的材料特性(如彈性模量、泊松比和密度等相關(guān)信息)和物理特性(主要是幾何參數(shù))。

邊界條件數(shù)據(jù)主要用以分析對(duì)象與外界的相互作用，主要包括位移約束和載荷條件。位移約束數(shù)據(jù)規(guī)定節(jié)點(diǎn)在自由度上所受到約束限制，載荷條件數(shù)據(jù)用于定義模型中節(jié)點(diǎn)載荷等力作用的位置、方向和大小。

同時(shí)，因?yàn)橛?jì)算程序中的數(shù)值計(jì)算大多為矩陣運(yùn)算，所以建立矩陣類class Matrix以及稀疏矩陣類class SparseMatrix，程序如下。

class Matrix

{

public:

Matrix();

Matrix(int,int);

Matrix(const Matrix &m);

Matrix(int,int,double);//預(yù)分配空間

virtual ～Matrix();//析構(gòu)函數(shù)

Matrix&operator=(const Matrix&);//矩陣的復(fù)制

int row() const;

int col() const;

private:

int rows_num,cols_num;

double **p;

void initialize();//初始化矩陣

};

矩陣類中部分函數(shù)定義如下。

double Point(int i,int j) const;

//求矩陣的逆矩陣

static Matrix inv(Matrix);

//矩陣加法

Matrix&operator+=(const Matrix&);

//矩陣減法

Matrix&operator-=(const Matrix&);

//矩陣乘法

Matrix&operator*=(const Matrix&);

//求K矩陣

Matrix&Kjuzhen(Matrix,Matrix);

class SparseMatrix

{

int nrLines,nrColumns;

int nrElements;//非零元素?cái)?shù)量

double* elements;//三元組的值存放

int* lines,*cols;//三元組的坐標(biāo)

public:

SparseMatrix();

SparseMatrix(int,int,double*,int*,int*);

SparseMatrix(int,int,int,double*,int*,int*);

SparseMatrix(const SparseMatrix&);

～SparseMatrix();

int getNumberOfLines();

int getNumberOfColumns();

int getElement();

}

稀疏矩陣中的部分函數(shù)如下。

friend SparseMatrix operator+(const

SparseMatrix&,const SparseMatrix&);

friend SparseMatrix operator-(const

SparseMatrix&,const SparseMatrix&);

friend SparseMatrix operator*(const

SparseMatrix&,const SparseMatrix&);

friend SparseMatrix operator*(const

SparseMatrix&,double)

這些函數(shù)實(shí)現(xiàn)計(jì)算程序全部的矩陣運(yùn)算，對(duì)矩陣的各種操作進(jìn)行封裝，能大大提高計(jì)算效率。

3.3.2 數(shù)值流程

彈塑性問題的改進(jìn)無(wú)單元Galerkin方法的數(shù)值計(jì)算流程如下。

(1)對(duì)結(jié)構(gòu)視為完全彈性，施加全部載荷進(jìn)行計(jì)算。

a.遍歷循環(huán)每個(gè)網(wǎng)格，對(duì)網(wǎng)格內(nèi)的高斯積分點(diǎn)進(jìn)行循環(huán)，判斷高斯積分點(diǎn)是否在域內(nèi)，是則進(jìn)行以下步驟，否則結(jié)束循環(huán)；確定高斯積分點(diǎn)影響域內(nèi)的節(jié)點(diǎn)，然后計(jì)算積分點(diǎn)處的形函數(shù)及其導(dǎo)數(shù)；用式(30)計(jì)算等效剛度矩陣，保存該積分點(diǎn)對(duì)矩陣及等效節(jié)點(diǎn)載荷的貢獻(xiàn)；

b.結(jié)束網(wǎng)格循環(huán)，得到剛度矩陣和等效載荷列向量Δ的第一項(xiàng)；

c.進(jìn)行邊界積分運(yùn)算，用類似方法計(jì)算等效載荷列向量Δ的第二項(xiàng)；

d.由第一項(xiàng)和第二項(xiàng)計(jì)算Δ；

e.擬合各節(jié)點(diǎn)的位移、應(yīng)力和應(yīng)變等信息。

(3)施加載荷增量Δ；

(4)按照式(25)計(jì)算剛度矩陣，要保證在彈性階段和塑性階段分別采用不同公式進(jìn)行計(jì)算；

(5)求解式(30)，得到所有節(jié)點(diǎn)的位移增量，進(jìn)而計(jì)算應(yīng)變?cè)隽亢蛻?yīng)力增量，并把位移增量、應(yīng)變?cè)隽亢蛻?yīng)力增量疊加到加載前的數(shù)據(jù)上；

(6)若加載到全部載荷，則跳出循環(huán)；

(7)輸出節(jié)點(diǎn)的位移。

4 數(shù)值算例

對(duì)二維和三維受集中力的懸臂梁采用彈塑性問題的改進(jìn)無(wú)單元Galerkin方法進(jìn)行彈塑性大變形數(shù)值分析。對(duì)懸臂梁建模，利用前處理程序模塊在求解域中布置節(jié)點(diǎn)，并將節(jié)點(diǎn)信息輸入到編寫的C++數(shù)值算例程序，計(jì)算獲取節(jié)點(diǎn)處的位移和應(yīng)力、應(yīng)變數(shù)據(jù)，并將計(jì)算結(jié)果與使用MATLAB軟件得到的數(shù)值解進(jìn)行比較，驗(yàn)證改進(jìn)無(wú)單元Galerkin方法以及程序的有效性。

4.1 二維算例

自由端受集中力載荷的二維懸臂梁示意見圖2。梁的幾何尺寸為=8 m，=1 m，載荷=1 N，泊松比=0.25，屈服極限=25 Pa。材料的彈性模量為=1.0×10Pa，按平面應(yīng)力問題求解。

圖2 受集中力的二維懸臂梁示意

采用前處理程序布置節(jié)點(diǎn)，結(jié)果見圖3。

圖3 二維懸臂梁節(jié)點(diǎn)分布

采用彈塑性問題的改進(jìn)無(wú)單元Galerkin方法進(jìn)行計(jì)算，將獲得的位移信息輸入后處理程序中，懸臂梁的變形見圖4，其位移量增大30倍。

圖4 二維懸臂梁變形

為驗(yàn)證程序的準(zhǔn)確性，將本文程序計(jì)算得到的數(shù)值解與MATLAB結(jié)果進(jìn)行對(duì)比，見圖5。本文二維計(jì)算程序能夠得到滿意的結(jié)果。

圖5 二維懸臂梁的節(jié)點(diǎn)位移計(jì)算結(jié)果對(duì)比

4.2 三維算例

自由端受集中力載荷的三維懸臂梁示意見圖6。梁的幾何尺寸為=8 m，=1 m，=1 m，其他材料特性與第4.1節(jié)相同，本文案例可采用線性強(qiáng)化彈塑性模型和MISES屈服準(zhǔn)則。

圖6 受集中力三維懸臂梁示意

采用前處理程序布置節(jié)點(diǎn)，結(jié)果見圖7。將計(jì)算得到的數(shù)據(jù)輸入后處理程序中，模型的變形結(jié)果見圖8和9，其位移量增大15倍。

圖7 三維懸臂梁節(jié)點(diǎn)分布

圖8 三維懸臂梁實(shí)體變形

圖9 三維懸臂梁節(jié)點(diǎn)變形

將本文程序計(jì)算得到的數(shù)值解與MATLAB結(jié)果對(duì)比，見圖10。本文軟件的三維計(jì)算程序能夠得到滿意的結(jié)果。

圖10 三維懸臂梁的節(jié)點(diǎn)位移計(jì)算結(jié)果對(duì)比

5 結(jié) 論

采用C++語(yǔ)言開發(fā)彈塑性問題的改進(jìn)無(wú)單元Galerkin方法計(jì)算軟件，實(shí)現(xiàn)快速建模和對(duì)離散點(diǎn)集進(jìn)行前處理后的數(shù)據(jù)采集，采用改進(jìn)無(wú)單元Galerkin方法的程序計(jì)算結(jié)果，將得到的節(jié)點(diǎn)信息和邊界條件等信息以文件形式輸入到后處理程序中，實(shí)現(xiàn)結(jié)果圖形可視化。通過(guò)數(shù)值算例分析，驗(yàn)證本文方法的有效性和優(yōu)越性。

在后續(xù)工作中，在前處理程序中將加入更多程序模塊，完善模型的邊界條件處理和組件功能等程序；計(jì)算程序模塊將加入彈性大變形和彈塑性大變形等問題的計(jì)算程序，以支持復(fù)雜問題的數(shù)值分析，為形成完善的無(wú)網(wǎng)格方法計(jì)算軟件奠定基礎(chǔ)。