第三章 線彈性靜力學(xué)問題

1、第三章 線彈性靜力學(xué)問題第一節(jié) 用向?qū)ppWizard解線彈性靜力學(xué)問題本例將求解一個(gè)簡支的深梁上表面受均布載荷作用下的平面應(yīng)變問題，求出梁的變形和應(yīng)力分布。圖3.1.1 受均布?jí)毫Φ暮喼盍涸搯栴}的偏微分方程表達(dá)式及定解條件如下：1）平衡方程 2）幾何方程 3）本構(gòu)方程 4）邊界條件第二類邊界條件具體操作步驟如下：設(shè)定工作路徑。點(diǎn)擊File菜單下的WorkDir選項(xiàng)，在彈出的對話框中填入工作路徑。如圖3.1.2所示圖3.1.2 設(shè)定路徑對話框點(diǎn)擊AppWizard菜單，彈出公式庫列表，系統(tǒng)提供了數(shù)學(xué)、固體力學(xué)、巖土、熱傳導(dǎo)、滲流、流體力學(xué)等領(lǐng)域的各種問題的公式。如圖3.1.3所示。點(diǎn)擊so

2、lid，在Project Name中填入“S”；點(diǎn)擊“Next”進(jìn)入下一步。圖3.1.3 選擇學(xué)科領(lǐng)域選擇坐標(biāo)系。點(diǎn)擊左邊的“2dxy”，選中二維直角坐標(biāo)系，如圖3.1.4所示。再點(diǎn)擊“Next”進(jìn)入下一步。圖3.1.4 選擇坐標(biāo)系選擇問題類型，即該問題是線彈性還是彈塑性。點(diǎn)擊“elastic”，選擇線彈性，如圖3.1.5所示。點(diǎn)擊“Next”進(jìn)入下一步。圖3.1.5 選擇問題類型選擇方程類型，即方程是平面應(yīng)力問題還是平面應(yīng)變問題，是靜態(tài)問題還是動(dòng)態(tài)問題；對于動(dòng)態(tài)問題還可選擇Newmark格式求解。點(diǎn)擊“Static_Plane_Strain_Problem”，選擇選擇靜態(tài)平面應(yīng)力問題的方程，

3、如圖3.1.6所示。點(diǎn)擊“Next”進(jìn)入下一步。圖3.1.6 選擇方程類型選擇單元類型。如圖3.1.7所示。點(diǎn)中“q4”，再點(diǎn)擊中間的“Add”按鈕，加入到右邊方框表示采用四節(jié)點(diǎn)矩形單元。點(diǎn)擊“Next”進(jìn)入下一步。圖3.1.7 選擇單元類型選擇邊界單元。對于本算例，因?yàn)槭菓?yīng)力邊界條件，所以需要邊界單元。點(diǎn)中“l(fā)l2”，再點(diǎn)擊中間的“Add”按鈕，加入到右邊方框表示采用兩節(jié)點(diǎn)線單元。如圖3.1.8所示。點(diǎn)擊“Next”進(jìn)入下一步。圖3.1.8 選擇邊界單元選擇求解器。系統(tǒng)提供有：對稱求解器（sin）、非對稱求解器（nin）、Gauss-Seidel迭代求解器（gs）、逐次超松弛迭代求解器（so

4、r）、共軛梯度法求解器（cgm）、不完全LU分解求解器（ilu）等。點(diǎn)擊“sin”，選擇對稱求解器。如圖3.1.9所示。點(diǎn)擊“Next”進(jìn)入下一步。圖3.1.9 選擇求解器選擇存儲(chǔ)方式。對于有限元程序自動(dòng)生成系統(tǒng)網(wǎng)絡(luò)版（IFEPG）僅提供用外存的存儲(chǔ)方式。對FEPG3.0提供外存和內(nèi)存兩種存儲(chǔ)方式。選擇計(jì)算數(shù)據(jù)。該數(shù)據(jù)給出求解區(qū)域的有限元數(shù)據(jù)，如圖3.1.10所示。用戶應(yīng)注意：點(diǎn)擊左邊文件列表選擇數(shù)據(jù)時(shí)，看一下右邊方框中說明的方程類型和單元類型與你前面選擇的是否一致，不一致系統(tǒng)將給出提示，一致時(shí)才可以正確運(yùn)行。點(diǎn)擊“sdispe_q4”，點(diǎn)擊“Next”進(jìn)入下一步。圖3.1.10 選擇有限元計(jì)

5、算數(shù)據(jù)自動(dòng)生成全部有限元程序。點(diǎn)擊“Run”即可生成全部有限元程序。如圖3.1.11和圖3.1.12所示。圖3.1.11 點(diǎn)擊“Run”自動(dòng)生成圖3.1.12 全部有限元程序計(jì)算結(jié)果插入res01.bmp圖3.1.13 算例的有限元計(jì)算結(jié)果第二節(jié) 改變施加的載荷 生成全部有限元計(jì)算程序之后，在FEPG界面的左邊是WorkSpace列表，它將所有生成的文件以及生成這些文件所需的原始文件都列在其中。FEPG施加荷載的方式有三種，一是分布荷載，作為材料屬性給出，二是面力，通過邊界單元材料屬性給出，三是節(jié)點(diǎn)力，通過節(jié)點(diǎn)給出。改變施加的載荷需修改前處理中的DAT文件，以下分別按三種加載分別敘述具體如何修

6、改。a. 改變分布力的大小在WorkSpace列表中找到“preprocessor”文件夾，在修改文件DAT文件之前先看一看MTI文件，該文件定義了有限元計(jì)算中需要的所有數(shù)據(jù)表格。單擊“preprocessor”文件夾中的s.mti文件，在界面的右邊則顯示出該文件。找到mate段如下maten pe pv fx fyi4 E10.4 E10.4 E10.4 E10.4mate 定義為材料參數(shù)表格名；n表示材料類型號(hào)；pe 表示彈性模量；pv 表示poisson比；fx 表示x方向的分布力；fy 表示y方向的分布力；在了解了mate定義的材料參數(shù)意義后，再看DAT文件。單擊“preprocess

7、or”文件夾中的s.dat文件編輯該文件，找到mate信息段內(nèi)容如下mate1;1.0e10;0.25;0.;0.;對應(yīng)與s.mti文件，我們可以知道1 表示第一種材料，如果有多種材料可在下一行寫“2”及相應(yīng)的材料參數(shù)；1. 0e10 表示材料的彈性模量的值；0.25 表示材料的泊松比的值；第一個(gè)0. 表示x方向分布力的大??；第二個(gè)0. 表示y方向分布力的大??；弄清各個(gè)數(shù)值的意義我們就可以對它作相應(yīng)的修改了，比如假設(shè)在x方向有1.0e3大小的力則mate信息段改成如下形式mate1;1.0e10;0.25;1.0e3;0.;下面看一下改變了載后的計(jì)算結(jié)果。首先在s.dat上單機(jī)鼠標(biāo)右鍵彈出如下

8、命令mti s dat該mti的兩個(gè)輸入?yún)?shù)分別是s.mti 和 s.dat兩個(gè)文件，輸出是前處理程序prg.for和前處理顯示程序prt.for。運(yùn)行上述命令重新生成prg.for和prt.for。然后運(yùn)行run_batch_file文件夾中的s.bat重新的計(jì)算結(jié)果如下插入res02.bmpb. 改變面力同樣，我們先看一下s.mti文件。在s.mti文件中找到bf段如下bfn fx fyi4 E10.4 E10.4bf 定義為面力的參數(shù)表格名n 表示面力類型fx 面的切向方向力的分量的值fy 面的法向方向力的分量的值（參考局部標(biāo)架說明）再看s.dat中對應(yīng)的行，找到bf信息段如下bf1;0

9、.0;-1.0e6;1 表示第一種類型的面力0.0 表示切向方向上的力為零-1.0e6 表示法向方向上的力為1.0e6，方向和法向相反。下面給出兩種不同的面力，一是將-1.0e6改為1.0e6，這表示算例中的梁是受向上的拉力，同上面所述保存s.dat文件并用mti命令重新生成prg.for和prt.for，再運(yùn)行s.bat重新進(jìn)行計(jì)算，得計(jì)算結(jié)果如下插入res03.bmp二是將bf信息段改為如下形式bf1;1.0e6;0.0;即在切向方向加上剪力，方向方向不加力。同樣重新生成prg.for并進(jìn)行計(jì)算。計(jì)算結(jié)果如下插入res04.bmpc. 施加節(jié)點(diǎn)力節(jié)點(diǎn)力的修改同樣需在前處理中給出，讓我們看一

10、看“preprocessor”中的s.mti文件和s.dat文件。先看s.mti中id段和disp段id fn idu idvi6 i6 i6id 表示該段定義約束信息n 表示節(jié)點(diǎn)編號(hào)idu x方向的約束信息存入idu數(shù)組idv y方向的約束信息存入idv數(shù)組disp fn u vi6 e10.4 e10.4disp 表示該段定義給定位移或荷載信息n 表示節(jié)點(diǎn)編號(hào)u x方向的給定位移或荷載信息存入u數(shù)組v y方向的給定位移或荷載信息存入v數(shù)組再看s.dat文件中id信息段和disp信息段id&nn=(nx+1)*(ny+1); i;1;1;i=1,nn1; -1;-1;nx+1; 1;

11、-1;&n=nn;&nn=(nx+1)*(ny+1) 定義總節(jié)點(diǎn)數(shù)i;1;1;I=1,nn 以配對的方括號(hào)表示循環(huán)，即表示先把所有節(jié)點(diǎn)的約束賦值為1。FEPG系統(tǒng)規(guī)定1表示節(jié)點(diǎn)自由，-1表示節(jié)點(diǎn)被約束。1; -1;-1; 表示1節(jié)點(diǎn)在x和y兩個(gè)方向都被約束住。nx+1;1;-1; 表示nx+1節(jié)點(diǎn)僅在y方向被約束住，x方向是自由的。&n=nn 表示把最大節(jié)點(diǎn)后賦值給n，用戶應(yīng)注意加上此行，因?yàn)檫@將使生成的prg.for中把約束信息寫入文件的循環(huán)的上限為最大節(jié)點(diǎn)號(hào)，也就是保證將所有的節(jié)點(diǎn)約束信息寫入文件。再看disp信息段dispi;0.0;0.0;i=1,(nx+1)*

12、(ny+1)disp 表示該段定義給定位移和荷載信息i;0.0;0.0;i=1,(nx+1)*(ny+1) 以配對的方括號(hào)表示循環(huán)，將所有節(jié)點(diǎn)在x和y兩個(gè)方向的位移或荷載信息都賦零。FEPG系統(tǒng)作如下判斷：如果某節(jié)點(diǎn)的某個(gè)方向的約束信息是-1，即該節(jié)點(diǎn)在該方向被約束住，則在disp中給出的相應(yīng)的值表示給定位移；如果某節(jié)點(diǎn)的某個(gè)方向的約束信息是1，即該節(jié)點(diǎn)在該方向是自由的，則在disp中給出的相應(yīng)的值表示給定荷載。為了更清楚地看到施加節(jié)點(diǎn)荷載的結(jié)果，先將分布力和面力各方向上的荷載值置零，修改后的mate信息段和bf信息段如下mate1;1.0e10;0.25;0.0;0.;bf1;0.0;0.0

13、;在右上角節(jié)點(diǎn)沿水平方向向右施加大小為1.0e4的集中力，修改后的disp信息段如下dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*(ny+1);1.0e4;0.0;用mti命令重新生成前處理程序再計(jì)算后的如下結(jié)果插入res05.bmp在右上角節(jié)點(diǎn)沿垂直方向向上施加大小為1.0e4的集中力，修改后的disp信息段如下dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*(ny+1); 0.0;1.0e4;用mti命令重新生成前處理程序再計(jì)算后的如下結(jié)果插入res06.bmp在左上角節(jié)點(diǎn)沿水平方向向左施加大小為1.0e4的集中力，修改后的disp

14、信息段如下dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*ny+1; -1.0e4; 0.0;&n=(nx+1)*(ny+1);大家應(yīng)注意此修改增加了一行賦值，即把總節(jié)點(diǎn)數(shù)賦給n，理由同在id中闡述的理由一樣。為何上兩例沒有給這樣的賦值呢？因?yàn)樯蟽衫际墙o的最大節(jié)點(diǎn)號(hào)，所以n自然取得總節(jié)點(diǎn)數(shù)，否則也應(yīng)有同樣的賦值行。荷載正負(fù)號(hào)的規(guī)定是和總體坐標(biāo)一致為正，反之為負(fù)。用mti命令重新生成前處理程序再計(jì)算后的如下結(jié)果插入res07.bmp, u1放大50000被的變形圖第三節(jié) 改變約束方式上面我們已經(jīng)提到id信息段給出約束信息，下面我們就舉幾個(gè)例子說明如何修改約束

15、信息假設(shè)物體左端固定，右端懸空，且在右上角受一大小為1.0e4向下的集中力，如下圖所示。對s.dat文件中id信息段做如下修改id&nn=(nx+1)*(ny+1);i;1;1;i=1,nn(i-1)*(nx+1)+1;-1;-1;i=1,ny+1&n=nn;相應(yīng)的節(jié)點(diǎn)荷載也作響應(yīng)修改，修改后內(nèi)容如下dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*(ny+1);0.0;-1.0e4;用mti命令重新生成前處理程序再計(jì)算后的如下結(jié)果插入res08.bmp，變形放大5000倍的結(jié)果假設(shè)物體左端固定，右端有一豎向支撐，且在右上角受一大小為1.0e4向下的集

16、中力，如下圖所示。在上例的基礎(chǔ)上只要增加右下角節(jié)點(diǎn)在y方向上的約束即可，修改后的id信息段如下id&nn=(nx+1)*(ny+1);i;1;1;i=1,nn(i-1)*(nx+1)+1;-1;-1;i=1,ny+1nx+1;1;-1;&n=nn;用mti命令重新生成前處理程序再計(jì)算后的如下結(jié)果插入res09.bmp，變形放大50000倍的結(jié)果最后一例假設(shè)物體兩端固定，且在中間位置受一大小為1.0e4向下的集中力，如下圖所示。修改id信息段如下id&nn=(nx+1)*(ny+1);i;1;1;i=1,nn(i-1)*(nx+1)+1;-1;-1;i=1,ny+1i*(

17、nx+1);-1;-1;i=1,ny+1&n=nn;disp信息段也作如下修改dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*ny+(nx+2)/2;0.0;-1.0e4;&n=(nx+1)*(ny+1);用mti命令重新生成前處理程序再計(jì)算后的如下結(jié)果插入res10.bmp，變形放大50000倍的結(jié)果第四節(jié) 改變材料參數(shù)前面在修改分布力是已經(jīng)提到mate信息段，并解釋了各參數(shù)的意義，其中前兩個(gè)參數(shù)是線彈性材料的彈性模量個(gè)泊松比，下面就修改材料參數(shù)，看計(jì)算結(jié)果如何變化。在兩段固定，中間受以集中力的情況下，看不同的材料的變形情況。下面分別給出在彈性模量

18、E分別為1.0e9，5.0e9和1.0e10時(shí)mate信息段修改后的內(nèi)容mate1;1.0e9;0.25;0.0;0.;mate1;5.0e9;0.25;0.0;0.;mate1;1.0e10;0.25;0.0;0.;每次修改后保存s.dat文件，用mti命令重新生成前處理程序再計(jì)算后的得如下計(jì)算結(jié)果的比較插入res11.bmp，變形放大5000倍的結(jié)果第五節(jié) 兩種不同的材料III下面考慮兩種不同材料時(shí)，如何用FEPG求解，假設(shè)求解下列兩層材料深梁在集中力作用下的變形，如下圖所示。如圖中的兩層材料彈性模量，泊松比，厚度相同，右上角受垂直向下大小為1.0e4大小的力，分析在此情景下深梁的變形。為

19、用FEPG解算該問題，需在前處理中作相應(yīng)的修改。看一看“preprocessor”文件夾中s.mti文件的q4段，內(nèi)容如下q4n nod1 nod2 nod3 nod4 matei6 i6 i6 i6 i6 i6q4 表示四節(jié)點(diǎn)矩形單元n 表示單元編號(hào)nod1 表示單元第一節(jié)點(diǎn)編號(hào)存入的數(shù)組名nod2 表示單元第二節(jié)點(diǎn)編號(hào)存入的數(shù)組名nod3 表示單元第三節(jié)點(diǎn)編號(hào)存入的數(shù)組名nod4 表示單元第四節(jié)點(diǎn)編號(hào)存入的數(shù)組名mate 表示該單元對應(yīng)的材料號(hào)對于本算例需修改s.dat文件實(shí)現(xiàn)分層劃分單元，上層mate的值給1，下層mate的值給2，mate對應(yīng)的信息段分別給出兩中材料的材料參數(shù)，修改后的

20、具體內(nèi)容如下。對q4信息段作如下修改q4nx*(j-1)+i;(j-1)*(nx+1)+i;(j- 1)*(nx+1)+i+1;j*(nx+1)+i+1;j*(nx+1)+i;1;i=1,nx;j=1,ny/2nx*(j-1)+i;(j-1)*(nx+1)+i;(j- 1)*(nx+1)+i+1;j*(nx+1)+i+1;j*(nx+1)+i;2;i=1,nx;j=ny/2+1,nymate信息段也需相應(yīng)修改mate1;1.0e10;0.25;0.0;0.;2;2.0e10;0.30;0.0;0.;用mti命令重新生成前處理程序再計(jì)算后得如下結(jié)果insert res12.bmp and res

21、13.bmp 變形放大5000倍的結(jié)果第六節(jié) 線彈性力學(xué)方程的PDE文件FEPG通過用戶填寫的PDE文件生成單元子程序。在了解如何填寫PDE文件之前先回顧一下線彈性力學(xué)所對應(yīng)的偏微分方程。線彈性力學(xué)偏微分方程包括平衡方程、幾何方程和本構(gòu)方程。在二維直接坐標(biāo)系下，線彈性平面應(yīng)變問題的方程具體可寫成如下形式平衡方程： （1）幾何方程 （2）本構(gòu)方程 （3）對(1)第一個(gè)方程乘，第二個(gè)方程乘，將兩式對體積積分并相加得 （4）對(4)式分布積分可得 （5）在考慮下列Cauchy公式 （6）其中和分別表示應(yīng)力邊界上力矢量在x和y方向的分量。將Cauchy公式和幾何方程（2）式代入（5）式得 （7）將本構(gòu)

22、方程代入（7）式得最終的虛功方程的弱形式 （8）其中 為一常系數(shù)有了(8)式我們就可以開始填寫PDE文件，先考慮（8）式中的體積積分項(xiàng)，由此構(gòu)造體單元。至于(8)中的邊界積分項(xiàng)我們在下一節(jié)給出詳細(xì)說明。按照（8）式給出如下的PDE文件，參閱“sdispe”文件夾中的sdispe.pde文件右邊給出相應(yīng)的注解disp u,v, 給出未知函數(shù)名U，Vcoor x,y, 給出總體坐標(biāo)系下的坐標(biāo)變量名X，Yfunc = funa,funb,func, 給出需要用到的函數(shù)名Shap %1 %2 給出單元形狀類型符和節(jié)點(diǎn)個(gè)數(shù)gaus %3 給出每個(gè)方向上積分點(diǎn)的個(gè)數(shù)或單元形狀類型符mass %1 給出單元

23、形狀類型符load = fu fv 給出載荷向量，對應(yīng)于$c6 pe = prmt(1) 從PRMT數(shù)組中取出彈性模量，$c6表示插入FORTRAN源程序$c6 pv = prmt(2) 從PRMT數(shù)組中取出泊松比$c6 fu = prmt(3) 從PRMT數(shù)組中取出X方向體力$c6 fv = prmt(4) 從PRMT數(shù)組中取出Y方向體力$c6 fact = pe/(1.+pv)/(1.-2.*pv) 由彈模和泊松比組合成系數(shù)FACT Func 給出一些函數(shù)的定義funa=+u/x 定義FUNA為，u/x表示 funb=+v/y 定義FUNB為，v/y表示 func=+u/y+v/x 定義F

24、UNC為stif 給出剛度矩陣dist =+funa;funa*fact*(1.-pv) funa;funa的第一個(gè)FUNA+funa;funb*fact*(pv) 表示，第二個(gè)FUNA表示+funb;funa*fact*(pv) ；表示積分。+funb;funb*fact*(1.-pv) +func;func*fact*(0.5-pv)load =+u*fu+v*fv 關(guān)于PDE文件還需作如下說明a. 單元形狀類型符，q表示矩形單元；b. guas信息段后給的參數(shù)是數(shù)字表示每個(gè)方向上積分點(diǎn)的個(gè)數(shù)，如后跟單元形狀類型符則表示采用節(jié)點(diǎn)積分；c. prmt( )數(shù)組保存了我們前面所說的前處理數(shù)據(jù)中

25、的材料參數(shù)，在PDE文件里通過賦值取出來。在了解了PDE文件如何得來以及PDE文件各內(nèi)容的意義后，讓我們看一看如何由PDE文件生成單元子程序，又如何連接到單元計(jì)算的E程序中。將sdispe.pde文件中的shap、gaus、mass信息段進(jìn)行以下修改shap q 4gaus qmass q其它內(nèi)容不變。上述修改表示將生成四節(jié)點(diǎn)矩形單元子程序，體積分采用節(jié)點(diǎn)積分。保存sdispe.pde文件，單擊FEPG界面主菜單中的“Run”，彈出FEPG命令對話框，在輸入欄輸入pdeges sdispe seuq4其中pdeges為系統(tǒng)命令，其功能是由PDE文件生成GES文件和單元子程序，它的第一個(gè)參數(shù)為P

26、DE文件主名，第二個(gè)參數(shù)為將要生成的GES文件也是單元子程序文件的主名。第二個(gè)參數(shù)也可不寫，此時(shí)表示生成的GES文件或單元子程序文件和PDE文件有相同的主名。另PDE文件和GES文件的關(guān)系是：將PDE文件中的shap段和gaus段用形函數(shù)庫和高斯積分庫作替換即得GES文件。上述命令表示由sdispe.pde文件生成seuq4.ges和seuq4.for。輸入上述命令后回車，F(xiàn)EPG即可生成用戶所需的單元子程序。把新生成的單元子程序連接到單元計(jì)算的E程序中，使用FEPG的let命令，具體做法是在“sdispe”文件夾上單擊鼠標(biāo)右鍵彈出let命令，選擇并單擊該命令即可實(shí)現(xiàn)把單元子程序連接到單元計(jì)算

27、的E程序中。第七節(jié) 邊界單元在上面講解線彈性力學(xué)方程PDE文件時(shí)，我們得到了線彈性力學(xué)虛功方程的弱形式。上面所說的PDE文件的給出了虛功方程弱形式中體積積分的描述，這里我們將考慮如何描述邊界積分項(xiàng)。FEPG是通過添加邊界單元來描述邊界積分的，邊界單元的作用相當(dāng)于加力的邊界條件或彈簧邊界，這就等效于偏微分方程的第二類或第三類邊界。下面我們看邊界單元的PDE文件。disp u,v, 給出未知函數(shù)名u，vcoor x, 給出單元局部坐標(biāo)系下的坐標(biāo)變量名Shap l 2 給出單元形狀類型符和節(jié)點(diǎn)個(gè)數(shù)，l 表示線單元， 2 表示兩節(jié)點(diǎn)單元gaus 2 給出每個(gè)方向上積分點(diǎn)的個(gè)數(shù)$c6 fx=prmt(1

28、) 從前處理中取出切向表面力， $c6表示插入FORTRAN源程序$c6 fy=prmt(2) 從前處理中取出法向表面力Stif 給出剛度矩陣dist=u;u*0.0 此算例為力的邊界，所以這里剛度矩陣等于零load=+u*fx+v*fy 給出載荷向量，對應(yīng)于end 結(jié)束符同樣使用pdeges命令由PDE文件生成相應(yīng)的單元子程序。假設(shè)上面的文件名為loadbnd.pde，單擊FEPG界面主菜單中的“Run”，彈出FEPG命令對話框，在輸入欄輸入pdeges loadbnd seull2按回車，系統(tǒng)即可由loadbnd.pde文件生成seull2.ges文件和邊界單元的單元計(jì)算的FORTRAN源

29、程序文件seull2.for。把新生成的單元子程序連接到單元計(jì)算的E程序中，同樣使用FEPG的let命令，具體做法是在“sdispe”文件夾上單擊鼠標(biāo)右鍵彈出let命令，選擇并單擊該命令即可實(shí)現(xiàn)把單元子程序連接到單元計(jì)算的E程序中。下面我們再更進(jìn)一步考慮有彈簧邊界的情景，即有等效于偏微分方程的第三類邊界條件的邊界。此時(shí)邊界力矢量和和邊界上的位移有關(guān)，假設(shè) （9）其中，分別為切向和方向的彈簧剛度。將和的表達(dá)式代入（8）式得有第三類邊界的弱虛功方程（10）不難看出（10）式比（8）多一項(xiàng)，此即由第三類邊界引起的剛度矩陣的修正。下面我們再給出有第三類邊界的邊界單元的PDE文件的填寫disp u,v,

30、 給出未知函數(shù)名u，vcoor x, 給出單元局部坐標(biāo)系下的坐標(biāo)變量名Shap l 2 給出單元形狀類型符和節(jié)點(diǎn)個(gè)數(shù)，l 表示線單元， 2 表示兩節(jié)點(diǎn)單元gaus 2 給出每個(gè)方向上積分點(diǎn)的個(gè)數(shù)$c6 fx=prmt(1) 從前處理中取出切向表面力， $c6表示插入FORTRAN源程序$c6 fy=prmt(2) 從前處理中取出法向表面力$c6 ax=prmt(3) 從前處理中取出彈簧切向剛度$c6 ay=prmt(4) 從前處理中取出彈簧法向剛度Stif 給出剛度矩陣dist=u;u*ax+v;v*ay 給出剛度矩陣的修正，對應(yīng)于load=+u*fx+v*fy 給出載荷向量，對應(yīng)于end 結(jié)

31、束符用FEPG的pdeges命令就可由此文件生成有彈簧邊界的邊界單元子程序，再用let命令將單元子程序連接到E程序。這里還有一點(diǎn)需要說明的是需在前處理中多定義兩個(gè)參數(shù)，具體做法是：a. 修改前處理中s.mti文件bf段，修改后的內(nèi)容如下bfn fx fy ax ayi4 E10.4 E10.4 E10.4 E10.4b. 修改前處理中s.dat文件bf段，修改后的內(nèi)容如下bf1;0.0;0.0;1.0e8;2.0e8; 后兩個(gè)給出彈簧切向剛度和法向剛度下圖顯示第5節(jié)中算例在加入彈簧邊界的結(jié)果插入res14.bmp 放大5000后的變形圖第八節(jié) 坐標(biāo)變換邊界單元定義的未知函數(shù)是局部坐標(biāo)下的未知函

32、數(shù)，而最后求解的整體坐標(biāo)下的未知函數(shù)，所以需把局部坐標(biāo)下的未知函數(shù)轉(zhuǎn)換到整體坐標(biāo)下再求解。假設(shè)在邊界單元中定義的局部坐標(biāo)下的位移為，整體坐標(biāo)下的位移是，則用下式將局部坐標(biāo)下的位移轉(zhuǎn)變到整體坐標(biāo)下的位移 （11）其中為坐標(biāo)轉(zhuǎn)換矩陣。在程序中是如何實(shí)現(xiàn)這種坐標(biāo)變換的呢？FEPG在用戶添好邊界單元的PDE文件并生成相應(yīng)的單元子程序之后，通過getglt命令即可自動(dòng)生成坐標(biāo)變換的glt文件。下面對命令的使用加以說明getglt命令的格式如下getglt ls gs dim T第一個(gè)參數(shù)ls表示單元子程序文件名(ges文件主名)第二個(gè)參數(shù)gs表示坐標(biāo)變換文件名（glt文件主名）第三個(gè)參數(shù)dim表示轉(zhuǎn)換到

33、的整體坐標(biāo)系的維數(shù)，添入數(shù)字第四個(gè)參數(shù)T表示作坐標(biāo)變換，不添該參數(shù)表示不作坐標(biāo)變換，此情況出現(xiàn)在未知函數(shù)為標(biāo)量的情況。對于第7節(jié)中的邊界單元可以用下列命令生成坐標(biāo)變換子程序getglt seull2 seugl2 2 T此命令生成seugl2.glt和seugl2.for文件，其中seugl2.glt文件內(nèi)容如下defigsub = seugl2lsub = seull2gdim = 2ldim = 1gvar gu1,gv1,gu2,gv2,lvar lu1,lv1,lu2,lv2,node = 2vartlu1 = T1 gu1 gv1 lu2 = T1 gu2 gv2 lv1 = T2

34、gu1 gv1 lv2 = T2 gu2 gv2然后使用FEPG的glt命令由glt文件生成相應(yīng)的坐標(biāo)變換子程序。具體做法是在glt文件上單擊鼠標(biāo)右鍵，彈出并使用glt命令，或者單擊主菜單“Run”彈出FEPG命令列表對話框，在輸入欄輸入glt seugl2按回車后即可生成相應(yīng)的坐標(biāo)變換子程序。同樣需要使用let命令將它連接到單元計(jì)算的E程序中。第九節(jié) 連接的接口IO文件在前面我們已多次提到用let命令把單元子程序連接到單元計(jì)算的E程序中，這里我們介紹在實(shí)現(xiàn)連接中一個(gè)很重要的文件IO文件。我們結(jié)合在第1節(jié)給出的算例中計(jì)算位移文件夾sdispe中的sdispe.io文件來說明該文件的意義。sdi

35、spe.io文件的內(nèi)容如下，右邊文字給出相應(yīng)的注解。y表示有材料性質(zhì)表格2 表示有限元程序采用兩種單元seuq4 給出第一種單元的子程序名，對應(yīng)于subroutine seuq4()seuq4 給出第一種單元的模塊名，對應(yīng)于seuq4.objseugl2 給出第二種單元的子程序名，對應(yīng)于subroutine seugl2()seugl2,seull2 給出第二種單元的模塊名，對應(yīng)于seugl2.obj和seull2.obj上面已經(jīng)說過，第二種單元之所以連接兩個(gè)模塊是因?yàn)檫吔鐔卧窃诰植孔鴺?biāo)下定義的，需用一個(gè)模塊實(shí)現(xiàn)向整體坐標(biāo)的轉(zhuǎn)換。此時(shí)我們看一看E程序esdispe.for的調(diào)用情況，esdi

36、spe.for含有如下行 goto (1,2), ityp1 call seuq4(r,coef,prmt,Es,Em,Ec,Ef,ne) goto 32 call seugl2(r,coef,prmt,Es,Em,Ec,Ef,ne) goto 33 continue可以看出該E程序中調(diào)用了seuq4()和seugl2()子程序，而在含有seugl2()子程序中又有 call seull2(y,coefr,prmt,els,elm,eld,ell,num)即對seull2()的調(diào)用，為實(shí)現(xiàn)這樣的調(diào)用，在生成可執(zhí)行文件時(shí)需連接seuq4.obj、seugl2.obj和seull2.obj三個(gè)目標(biāo)

37、文件。下面由第六節(jié)中的PDE文件sdispe.pde，用下列命令生成另一個(gè)單元子程序。單擊Run主菜單，在命令輸入欄輸入pdeges sdispe seuq4my這樣由sdispe.pde生成了seuq4my.ges和seuq4my.for，并將seuq4my.for由FORTRAN編譯器編譯成sueq4my.obj。在sdispe文件夾上單擊鼠標(biāo)右鍵并選擇Add files to folder命令，加入seuq4my.ges和seuq4my.for可以看到其內(nèi)容和sueq4一樣，只是子程序名和for文件主名不一樣。下面修改sdispe.io，連接新生成的單元子程序。修改后的sdispe.io

38、文件如下y2seuq4myseuq4myseugl2seugl2,seull2在sdispe.nfe文件上單擊鼠標(biāo)右鍵，彈出并使用gnfepg命令重新生成E程序和U程序。生成后再次查看E程序esdispe.for中單元子程序的調(diào)用發(fā)生變化。 goto (1,2), ityp1 call seuq4my(r,coef,prmt,Es,Em,Ec,Ef,ne) goto 32 call seugl2(r,coef,prmt,Es,Em,Ec,Ef,ne) goto 33 continue所以用戶可以通過修改io文件實(shí)現(xiàn)對子程序的靈活調(diào)用。第十節(jié) 計(jì)算應(yīng)力的PDE文件應(yīng)力計(jì)算的單元子程序同樣也是由相

39、應(yīng)的PDE文件生成的。在講應(yīng)力計(jì)算單元的PDE文件如何填寫之前先介紹一下應(yīng)力的計(jì)算方法。大家知道，由位移為基本未知量得到的位移解在全域是連續(xù)的，應(yīng)力和應(yīng)變在單元內(nèi)部是連續(xù)的，而單元間是不連續(xù)的，即在單元邊界上發(fā)生突跳。因此，對同一個(gè)節(jié)點(diǎn)，圍繞它的不同的單元計(jì)算得到的應(yīng)力的值是不同的。FEPG采用最小二乘算法對由位移算得的應(yīng)力進(jìn)行局部磨平。構(gòu)造下列泛函 （12）其中為磨平后的應(yīng)力，表示由已知位移計(jì)算得到的單元應(yīng)力。 （13）其中表示彈性矩陣，表示已知位移由幾何方程求得的應(yīng)變。于是，用最小二乘算法對應(yīng)力進(jìn)行局部磨平等價(jià)于求（12）定義的泛函的極值。即求使 （14）而(14)式又等價(jià)于，即有 （15

40、）將（15）式改寫為 (16)至此我們可以填寫PDE文件對(16)進(jìn)行求解，計(jì)算應(yīng)力的PDE文件如下（右邊是注解）coef u v 位移作為已知函數(shù)傳入應(yīng)力計(jì)算單元coor x y 定義坐標(biāo)變量名disp sa,sb,sc, 定義待求的應(yīng)力的三個(gè)分量名shap %1 %2 給出單元類型符和節(jié)點(diǎn)數(shù)gaus %3 給出每個(gè)方向上積分點(diǎn)的個(gè)數(shù)或單元形狀類型符mass %1 給出單元形狀類型符，若系數(shù)不為1，需在第二個(gè)參數(shù)給出load = fsa fsb fsc 給出應(yīng)力函數(shù)名$c6 pe = prmt(1) 給出彈性模量，以下參數(shù)意義同位移計(jì)算的PDE文件$c6 pv = prmt(2)$c6 fx

41、 = prmt(3)$c6 fy = prmt(4)$c6 fact = pe/(1.+pv)/(1.-2.*pv)stif$cv funa=+u/x 計(jì)算x方向正應(yīng)變，*/*表示對已知函數(shù)的導(dǎo)數(shù)$cv funb=+v/y 計(jì)算y方向正應(yīng)變$c6 fsa=+funa*(1.-pv)+funb*(pv) $c6 fsa=fsa*fact 計(jì)算x方向的正應(yīng)力$c6 fsb=+funa*(pv)+funb*(1.-pv)$c6 fsb=fsb*fact 計(jì)算y方向的正應(yīng)力$cv func=+u/y+v/x 計(jì)算剪應(yīng)變$c6 fsc=func*(0.5-pv)$c6 fsc=fsc*fact 計(jì)算剪應(yīng)

42、力dist = sa;sa*0.0 系數(shù)矩陣為對角陣，在mass信息段給出，所以此處dist為零11計(jì)算主應(yīng)力如何使用FEPG輸出用戶定義的物理量，例如主應(yīng)力、按材料強(qiáng)度理論的折算應(yīng)力等，下面我們通過一個(gè)平面應(yīng)力問題說明如何輸出并顯示用戶定義的物理量。先用AppWizard生成一個(gè)平面應(yīng)力的例子，然后再做些修改，加上主應(yīng)力的計(jì)算，輸出并顯示計(jì)算結(jié)果。具體的操作步驟是：a. 設(shè)置新路徑；單擊File主菜單，再單擊Work Dir子菜單，彈出設(shè)置路徑對話框輸入新路徑或用瀏覽選擇一個(gè)已有路徑；b. 進(jìn)入FEPG向?qū)В粏螕鬉ppWizard主菜單，進(jìn)入向?qū)υ捒?；c. 選擇固體力學(xué)；在對話框左邊選擇S

43、olid，在Project Name欄輸入工程名比如“S”，單擊Next進(jìn)入下一步；d. 單擊2dxy，選擇二維直角坐標(biāo)系，單擊Next進(jìn)入下一步；e. 單擊elastic，選擇線彈性，單擊Next進(jìn)入下一步；f. 單擊Static_Plane_Stress_Problem，選擇靜態(tài)平面應(yīng)力問題，單擊Next進(jìn)入下一步；g. 單擊q4，再單擊Add按鈕，選擇四節(jié)點(diǎn)矩形單元，單擊Next進(jìn)入下一步；h. 單擊l2，再單擊Add按鈕，選擇兩節(jié)點(diǎn)線單元，單擊Next進(jìn)入下一步；i. 單擊sin，選擇對稱用內(nèi)存的求解器，單擊Next進(jìn)入下一步；j. 單擊outcore，選擇存儲(chǔ)單元?jiǎng)偠龋瑔螕鬘ext進(jìn)

44、入下一步；k. 單擊sdispd.q4，選擇好算例數(shù)據(jù)，單擊Next進(jìn)入下一步；l. 單擊Run，生成全部有限元程序；m. 單擊File主菜單，單擊Open Workspace，單擊文件夾中的wsp文件，顯示生成的程序。至此，我們點(diǎn)鼠標(biāo)的形式生成了，計(jì)算靜態(tài)線彈性平面應(yīng)力問題的全部有限元程序。下面介紹如何修改以加入主應(yīng)力。先看一下平面應(yīng)力問題計(jì)算應(yīng)力的PDE文件，右邊文字給出說明。最好與平面應(yīng)變問題計(jì)算應(yīng)力的PDE對照閱讀以加深理解。coef u v 位移作為已知函數(shù)傳入應(yīng)力計(jì)算單元coor x y 定義坐標(biāo)變量名disp sa,sb,sc, 定義待求的應(yīng)力的三個(gè)分量名shap %1 %2 給

45、出單元類型符和節(jié)點(diǎn)數(shù)gaus %3 給出每個(gè)方向上積分點(diǎn)的個(gè)數(shù)或單元形狀類型符mass %1 給出單元形狀類型符，若系數(shù)不為1，需在第二個(gè)參數(shù)給出load = fsa fsb fsc 給出應(yīng)力函數(shù)名$c6 pe = prmt(1) 給出彈性模量$c6 pv = prmt(2) 給出泊松比$c6 fu = prmt(3) 給出x方向上體積力$c6 fv = prmt(4) 給出y方向上體積力$c6 fact = pe/(1.+pv)/(1.-pv) 給出彈性矩陣中的一個(gè)系數(shù)，只是為了寫剛度信息 時(shí)方便，注意與平面應(yīng)變的差別stif$cv funa=+u/x 計(jì)算x方向正應(yīng)變，*/*表示對已知函數(shù)的導(dǎo)數(shù)$cv funb=+v/y 計(jì)算y方向正應(yīng)變$c6 fsa=+fu