第一章偏微分方程的“弱”形式

1、基于pFEPG有限元方法序 言本教材的目的是為了幫助pFEPG系統(tǒng)的用戶更好地理解和掌握該系統(tǒng)，從有限元方法原理，從有限元軟件的程序結構和數(shù)據(jù)結構，從有限元的輸入數(shù)據(jù)形式，以及軟件設計方法等方面對系統(tǒng)進行具體闡述。眾所周知，有限元方法是結構力學家發(fā)明的，其基本原理是基于變分法與分片插值多項式兩項技術，但是有限元技術的推廣與研究的深入使得今天的有限元方法已突破了這兩項技術的范圍。例如，為了把有限元方法應用于非結構力學領域，如流體力學領域，人們不得不采用虛位移原理即“弱形式”取代變分原理，為了求解無界區(qū)域問題或者奇點問題，采用非多項式的形函數(shù)（基函數(shù)）將具有更好的精度，因此形函數(shù)不得不突破分片插值

2、多項式的限制。pFEPG系統(tǒng)是基于虛位移原理（即弱形式）而不是變分原理，要求用戶書寫弱形式的微分方程表達式，因此本教材的第一章通過對不同領域的微分方程表達式推導其弱形式，讓用戶對微分方程弱形式的推導有一般性的了解，但是一個好的弱形式可能是專家經(jīng)過多年的研究后取得的成果，并不是輕而易舉就可得到的。教材的第二章給出各種常用的分片插值多項式形函數(shù)表達式及其推導的一般性方法，但本系統(tǒng)允許用戶使用非多項式形函數(shù)。第三章給出本系統(tǒng)有限元程序所需的輸入數(shù)據(jù)的一般形式。第四章論述了pFEPG系統(tǒng)的有限元計算程序結構，并給出全部FORTRAN源程序及其詳細的說明。第一章 偏微分方程的“弱”形式-虛位移原理11.

3、1 偏微分方程的弱解形式1 問題的提出1 偏微分方程弱解的積分形式-虛位移原理21.2 穩(wěn)態(tài)問題弱解的積分形式3 分部積分公式3 二維穩(wěn)態(tài)熱傳導問題的“弱”形式5 三維線彈性小變形靜態(tài)問題的“弱”形式8 三維穩(wěn)態(tài)滲流問題的“弱”形式14 二維不可壓縮流體穩(wěn)態(tài)Navier_Stokes方程的“弱”形式17 三維靜電場問題的“弱”形式20 三維柱坐標靜電場問題的“弱”形式221.3 瞬態(tài)問題的“弱”形式24 三維線彈性小變形動態(tài)問題的“弱”形式24 三維瞬態(tài)熱傳導問題的“弱”形式28 二維不可壓縮流體瞬態(tài)Navier_Stokes方程的“弱”形式311.4 求解解梯度的最小二乘法34 已知位移求應

4、力34 已知溫度求熱流密度36 已知電勢求電場強度37第二章 分片多項式的形函數(shù)392.1 插值函數(shù)與單元類型39 一維Lagrange單元39 二維單元41 三維單元492.2 等參單元55 導數(shù)之間的變換582.3 數(shù)值積分62 高斯積分63 節(jié)點積分65第三章 有限元輸入數(shù)據(jù)形式683.1 pFEPG系統(tǒng)的有限元輸入數(shù)據(jù)組成簡述68 輸入數(shù)據(jù)形式68 輸入數(shù)據(jù)框圖68 表格文件的讀寫格式693.2 單場問題的有限元輸入數(shù)據(jù)70 坐標數(shù)據(jù)表格70 節(jié)點規(guī)格數(shù)表格70 指定節(jié)點位移和節(jié)點荷載信息表格71 初始值表格71 單元信息數(shù)據(jù)723.3 有限元輸入數(shù)據(jù)的顯示和查詢733.4 pFEPG

5、系統(tǒng)的PRE文件74 線性的、與時間無關的問題74 非線性、依賴時間問題803.5 多場問題的有限元輸入數(shù)據(jù)87 場的命名約定88 多場問題舉例說明88第四章 有限元方法的源程序974.1 有限元程序結構與元件化程序設計方法97 程序結構97 元件化程序設計方法98 線性穩(wěn)態(tài)有限元問題99線性動態(tài)有限元問題100 非線性穩(wěn)態(tài)有限元問題101 非線性動態(tài)有限元問題1024.2 元件程序的結構104 元件程序104功能104 命令行參數(shù)說明105 參數(shù)及數(shù)組說明106 源程序106 Fortran源程序113 BFT元件程序114 功能114 命令行參數(shù)說明115 參數(shù)及數(shù)組說明115 源程序11

6、6 Fortran源程序120 E元件程序122 功能123 命令行參數(shù)說明124 參數(shù)及數(shù)組說明124 源程序125 Fortran源程序137 SOLV求解器139 功能139 命令行參數(shù)139 源程序139 Fortran源程序164 U元件程序167 功能167 命令行參數(shù)167參數(shù)及數(shù)組說明168 源程序168 Fortran源程序172第一章 偏微分方程的“弱”形式虛位移原理第一章 偏微分方程的“弱”形式虛位移原理1.1 偏微分方程的弱解形式 問題的提出 工程或物理學中的許多問題，通常是以未知場函數(shù)應滿足的偏微分方程和邊界條件的形式提出來的，可以一般地表示為未知函數(shù)應滿足偏微分方程

7、組 ()域可以是體積域、面積域等，如圖所示。同時未知函數(shù)還應滿足邊界條件 ()是域的邊界。 域 y 0 x圖要求解的未知函數(shù)可以是標量函數(shù)場(例如溫度)，也可以是幾個變量組成的向量函數(shù)場(例如位移、應變、應力等)。是表示對于獨立變量(例如空間坐標、時間坐標等)的微分算子。偏微分方程數(shù)應和未知場函數(shù)的數(shù)目相對應，因此，上述偏微分方程可以是單個的方程，也可以是一組方程。所以在式( 對于工程或物理學中遇到的偏微分方程一般是沒有理論解的，即未知函數(shù)沒有解析表達式。而工程上又需要了解這些未知函數(shù)，所以我們一般用數(shù)值的方法來求解這些偏微分方程。有限元方法就是一種數(shù)值求解偏微分方程的方法，它實際上求解的是偏

8、微分方程的弱解積分形式，所以我們需要先將偏微分方程變成其弱解積分形式，才能使用有限元方法。 偏微分方程弱解的積分形式-虛位移原理由于偏微分方程組()在域中每一點為零，因此就有 ()其中 ()是向量函數(shù)，我們稱為試驗函數(shù)或虛位移函數(shù)，它是一組和偏微分方程個數(shù)相等的任意函數(shù)。 我們稱式()對于任意的在域內某些點或部分子域中不滿足，即出現(xiàn)，馬上可以找到適當?shù)暮瘮?shù) 在很多情況下可以對()式進行分部積分得到另一種形式 ()其中，是微分算子，它們中所包含的未知函數(shù)導數(shù)的階數(shù)較()式的微分算子低，這樣對函數(shù)的連續(xù)性要求是以提高的連續(xù)性要求為代價的，由于原來對 (在(115)式中)并無連續(xù)性要求，但是適當提高

9、對其連續(xù)性的要求并不困難，因為它們是可以選擇的已知函數(shù)。這種降低對函數(shù)“弱”1.2 穩(wěn)態(tài)問題弱解的積分形式 偏微分方程一般可分為未知函數(shù)與時間有關的瞬（動）態(tài)問題和未知函數(shù)與時間無關的穩(wěn)（靜）態(tài)問題。在本節(jié)中，我們主要討論穩(wěn)態(tài)問題弱解的積分形式，下節(jié)討論瞬態(tài)問題弱解的積分形式。對于本節(jié)第一小節(jié)給出分部積分的一些基本公式，以下各節(jié)將用到這些公式。 分部積分公式在本節(jié)中，我們給出一些經(jīng)常用的分部積分公式。設和是三維區(qū)域上的類型標量函數(shù)，我們有如下的分部積分公式： () () ()其中為邊界上外法線方向的單位向量的三個分量，可寫成：式中，是指坐標軸方向與單位向量之間夾角的余弦。設和是三維區(qū)域上的類型

10、向量函數(shù)，其相應三個分量分別為、和、。令和分別代表三維直角坐標系下的梯度算子和Laplace算子，即令和分別代表三維直角坐標系下的散度算子和旋度算子，即由上面分部積分公式及算子定義，可推出下列等式成立，這些等式在把偏微分方程轉變?yōu)槠浞e分“弱”形式的過程中非常有用。 () () () 其中，點表示向量的標量積，表示向量的向量積；表示域的表面的法向單位向量；表示標準的法向梯度算子： 二維穩(wěn)態(tài)熱傳導問題的“弱”形式本小節(jié)通過一個穩(wěn)態(tài)熱傳導問題，我們來看看將偏微分方程化為其弱解積分形式的一般過程。對于二維直角坐標系，穩(wěn)態(tài)熱傳導方程如下 邊界條件如下這里表示溫度，是熱傳導系數(shù)，和是邊界上溫度和熱流的給定

11、值；是熱源密度乘以材料密度；由方程，我們寫出相當于（）的形式如下： ()其中，是任意標量函數(shù)。對式（ ()利用方程的邊界條件，()可變?yōu)椋?()若在選擇未知函數(shù)時，已滿足強制邊界條件即在邊界上滿足，則可以通過適當選擇，使在邊界上而略去(1210)式中沿邊界積分項，使相應的“弱”形式取得更簡潔的表達式如下： ()這樣我們就得到了最終的偏微分方程“弱”形式。對于一般問題步驟如下：1 先將偏微分方程化為其積分形式；2 利用分部積分公式將其積分形式化為“弱”形式；3 利用邊界條件將“弱”形式化為更簡潔的表達式。特別指出對于強制邊界條件（即第一類邊界條件），這時未知函數(shù)在此類邊界上的值已確定，可以選取虛

12、位移函數(shù)在此類邊界上的值為0，則 “弱”形式可略去沿此類邊界上的邊界積分項。有了偏微分方程的“弱”形式，我們就可以書寫有限元自動生成系統(tǒng)（pFEPG）的PDE(或VDE)文件，然后通過pFEPG的單元子程序自動生成系統(tǒng)來生成單元子程序的FORTRAN77源程序了。在本章中，我們將給出例子的PDE或VDE文件，并做簡單的說明，關于pFEPG的單元子程序自動生成系統(tǒng)的詳細說明請參見pFEPG說明書中關于PDE或VDE文件填寫的有關章節(jié)。偏微分方程等效積分的“弱”形式中一般有兩類積分項（一種為體積分項，一種為邊界積分項），所以我們需要填寫兩個PDE文件，對于邊界積分項的PDE文件，在pFEPG系統(tǒng)中

13、，我們稱為FBC文件。關于上面的例子中，填寫的pFEPG系統(tǒng)文件如下：PDE文件（右邊文字是注釋）disp u,coor x,y,shap %1 %2gaus %3mate ek eq 4.4e-2;0.0;0.0; ek，eq分別對應弱形式中的 空一行stifdist =+u/x;u/x*ek+u/y;u/y*ek 對應弱形式體左端積分項 空一行l(wèi)oad =+u*eq 對應弱形式右端體積分項 空一行endFBC文件（右邊文字是注釋）disp ucoor xshap %1 %2gaus %3mate eb 3.0e1; eb對應弱形式中的 空一行stifdist=+u;u*0.0 對應弱形式左

14、端邊界積分項，例子中沒有，所以系數(shù)取0.0 空一行l(wèi)oad=+u*eb 對應弱形式右端邊界積分項 空一行end 三維線彈性小變形靜態(tài)問題的“弱”形式對于三維直角坐標問題，線彈性小變形靜態(tài)問題基本方程如下：1.彈性體V域內任一點沿x,y,z方向的平衡方程為其中為表示體內任一點的應力狀態(tài)的6個應力分量，為正應力，為剪應力，為單位體積的體積力在x,y,z方向的分量。2.在微小位移和微小變形的情況下，略去位移導數(shù)的高次冪，則應變和位移間的幾何方程為：其中為表示體內任一點的應變的6個應變分量，為正應變，為剪應變，為體內任一點的在x,y,z方向的三個位移分量。3.對于各向同性的線彈性材料，應力和應變間的本

15、構方程為：其中是彈性體材料的彈性模量，是泊桑比。本構方程的矩陣表示為其中稱為彈性矩陣。為應力張量，為應變張量4.彈性體的全部邊界為。一部分邊界上已知外力。稱為力的邊界條件，這部分邊界用表示；另一部分邊界上彈性體的位移已知，稱為幾何邊界條件或位移邊界條件，這部分邊界用表示。這兩部分邊界構成彈性體的全部邊界，即力的邊界條件為：其中，為彈性體在邊界上單位面積的內力。設邊界外法線為N，其方向余弦為，則邊界上彈性體的內力可由下式確定：力邊界條件表示為其中為已知面力向量位移邊界條件為：其中，為彈性體在邊界上已知的位移。位移邊界條件表示為其中，為位移向量為已知位移向量。運用虛位移原理求位移由平衡方程，得其中

16、，為三個方向的虛位移。由本章節(jié)中的分部積分公式將其化為弱形式，利用幾何方程可得上式表示為向量形式如下其中，點表示向量的標量積。為體積力向量將本構方程的矩陣表示代入上式，可得將力和位移邊界條件代入上式（注意，對于位移邊界條件，虛位移為0），可得由上面的弱形式，填寫的pFEPG系統(tǒng)文件如下：VDE文件（右邊文字是注釋）disp u v wcoor x y zfunc exx eyy ezz eyz exz exymate pe pv fu fv fw 1.0d10;0.3;0.0;0.0;0.0; pe,pv分別對應弱形式中的, fu,fv,fw分別對應弱形式中的shap %1 %2gaus %3

17、vect u u v wvect f fu fv fwvect ev exx eyy ezzvect eg eyz exz exyvect dg (0.5-pv) (0.5-pv) (0.5-pv)matr dv 3 3(1.-pv) pv pvpv (1.-pv) pvpv pv (1.-pv) 空一行$c6 fact = pe/(1.+pv)/(1.-2.*pv) 空一行funcexx=+u/x 空一行eyy=+v/y 空一行ezz=+w/z 空一行eyz=+v/z+w/y 空一行exz=+w/x+u/z 空一行exy=+u/y+v/x 空一行stifdist=+ev_i;ev_j*dv_

18、i_j*fact+eg_i;eg_i*dg_i*fact 對應弱形式左端體積分項 空一行l(wèi)oad=+u_i*f_i 對應弱形式右端體積分項 空一行endFBC文件（右邊文字是注釋）disp u,v,w,coor x,y,mate fx fy fz 0.0;0.0;1.0d2; fx,fy,fz分別對應弱形式中的shap %1 %2gaus %3 空一行stifdist=+u;u*0.0 對應弱形式左端邊界積分項，例子中沒有，所以系數(shù)取0.0 空一行l(wèi)oad=+u*fx+v*fy+w*fz 對應弱形式右端邊界積分項 空一行end 三維穩(wěn)態(tài)滲流問題的“弱”形式對于三維直角坐標問題，穩(wěn)態(tài)滲流問題基本

19、方程如下：1對于不可壓縮的流體域內，連續(xù)方程為：其中，為方向的流速分量，為內源。2若介質是各向異性的，根據(jù)廣義達西定律，在方向的流速分量分別為：其中，為壓力函數(shù)，為九個滲透系數(shù)，由于對稱性，獨立系數(shù)共有6個。其矩陣表示為：其中，為滲流矩陣為流速向量3. 流體的全部邊界為。一部分邊界上單位面積滲流量已知,即法向流速已知其中，分別為邊界表面外法線在方向的方向余弦。另一部分邊界上壓力已知在工程上，這兩部分邊界經(jīng)常構成了流體的全部邊界，即利用虛位移原理求壓力由上面的連續(xù)方程可得：其中，為壓力的虛位移。由本章節(jié)中的分部積分公式將其化為弱形式，可得上式表示為向量形式如下將流速向量代入上式，可得將邊界條件代

20、入上式（注意，對于已知壓力邊界條件，虛位移為0），可得由上面的弱形式，填寫的pFEPG系統(tǒng)文件如下：VDE文件（右邊文字是注釋,注意VDE文件中的未知量用u表示）disp u coor x,y,zfunc exx eyy ezzmate pkx pky pkz pkxy pkyz pkxz eq pkx,pky,pkz,pkxy,pkyz,pkxz分別對應弱形式中的,eq對應弱形式中的shap %1 %2gaus %3vect ev exx eyy ezzmatr ek 3 3pkx pkxy pkxzpkxy pky pkyzpkxz pkyz pkz 空一行 空一行funcexx=+u/x

21、 空一行eyy=+u/y 空一行ezz=+u/z 空一行stifdist=+ev_i;ev_j*ek_i_j 對應弱形式左端體積分項 空一行l(wèi)oad=+u*eq 對應弱形式右端體積分項 空一行endFBC文件（右邊文字是注釋）disp ucoor x,y,zmate ev ev對應弱形式中的shap %1 %2gaus %3 空一行stifdist=+u;u*0.0 對應弱形式左端邊界積分項，例子中沒有，所以系數(shù)取0.0 空一行l(wèi)oad=+u*ev 對應弱形式右端邊界積分項 空一行end 二維不可壓縮流體穩(wěn)態(tài)Navier_Stokes方程的“弱”形式對于二維直角坐標問題，不可壓縮流體穩(wěn)態(tài)Nav

22、ier_Stokes方程如下：1.對于流體域內，連續(xù)性方程為 其中，為的速度分量。2.動量方程為其中，為流體壓力，分別為流體域內方向上的單位外力；為流體密度常數(shù)，為粘性系數(shù)，為一常數(shù)。關于邊界條件，我們將在下面推導弱形式時提出，便易理解。利用虛位移原理求解首先，利用連續(xù)性方程可得：其中為的虛位移。利用動量方程可得：其中，分別為的虛位移。由本章節(jié)中的分部積分公式將其化為弱形式，可得其中，為流體的全部邊界，將上面的三個弱形式合并可得：下面提出邊界條件，在流體進口處，邊界上的速度向量已知，即此時，虛位移為0。在固壁上，粘附邊界上速度向量為0，即此時，虛位移也為0；在流體出口處()， 將邊界條件代入上

23、面的弱形式，可得由上面的弱形式，填寫的pFEPG系統(tǒng)文件如下：VDE文件（右邊文字是注釋）disp u,v,p coor x,ycoef un vn pn 導入上一迭代步的u,v,p,用于對非線性項的處理func divmate rou emu fx fy rou,emu,fx,fy分別對應弱形式中的shap %1 %2gaus %3vect x x yvect f fx fyvect u u vvect un un vn 空一行funcdiv=+u_i/x_i 空一行stifdist= +u_i/x_j;u_i/x_j*emu -div;p-p;div+u/x_j;u*un_j*rou+v/

24、x_j;v*un_j*rou 對應弱形式左端體積分項 空一行l(wèi)oad=+u_i*f_i*rou 對應弱形式右端體積分項 空一行end這個問題屬于非線性問題，需要線性化才能計算，在這里我們用的是簡單迭代法進行線性化的。 三維靜電場問題的“弱”形式對于三維直角坐標問題，對于介質域內靜電場問題方程如下：1. 位勢方程其中，為電勢，是介質的介電常數(shù)，是介質內的體電荷密度。2. 介質的全部邊界為。一部分邊界為Dirichlet邊界條件，即規(guī)定了導體面上的電勢= 另一部分邊界上是齊次Neumann邊界條件即電勢法向導數(shù)為零其中，分別為邊界表面外法線在方向的方向余弦。它適用于對稱面上。在工程上，這兩部分邊界

25、經(jīng)常構成了介質的全部邊界，即運用虛位移原理求電勢利用位勢方程可得其等效積分形式為：其中，為電勢的虛位移。由本章節(jié)中的分部積分公式將其化為弱形式，可得將邊界條件代入上式（注意，對于Dirichlet邊界條件，虛位移為0），可得由上面的弱形式，填寫的pFEPG系統(tǒng)文件如下：PDE文件（右邊文字是注釋,注意PDE文件中的未知量用u表示）disp u,coor x,y,zshap %1 %2gaus %3mate ep rou 1.0;0.0; ep,rou分別對應弱形式中的 空一行stifdist=+u/x;u/x*ep+u/y;u/y*ep+u/z;u/z*ep 對應弱形式左端體積分項 空一行l(wèi)o

26、ad=+u*rou 對應弱形式右端體積分項 空一行end 三維柱坐標靜電場問題的“弱”形式對于三維柱坐標系下，介質域內靜電場問題方程如下：1. 位勢方程其中，為電勢，是介質的介電常數(shù)，是介質內的體電荷密度。2. 介質的全部邊界為。一部分邊界為Dirichlet邊界條件，即規(guī)定了導體面上的電勢= 另一部分邊界上是齊次Neumann邊界條件即電勢法向導數(shù)為零其中，分別為邊界表面外法線在方向的方向余弦。它適用于對稱面上。在工程上，這兩部分邊界經(jīng)常構成了介質的全部邊界，即梯度算子為運用虛位移原理求電勢利用位勢方程可得其等效積分形式為：其中，為電勢的虛位移。由本章節(jié)中的分部積分公式將其化為弱形式，可得將

27、邊界條件代入上式（注意，對于Dirichlet邊界條件，虛位移為0），可得由上面的弱形式，填寫的pFEPG系統(tǒng)文件如下：PDE文件（右邊文字是注釋,注意PDE文件中的未知量用u表示）disp u,coor r,o,zshap %1 %2gaus %3mate ep rou 1.0;0.0; ep,rou分別對應弱形式中的$i$c6 r0=0.0$c6 do 1001 n=1,nnode$c6 r0=r0+coorr(1,n)$c0 1001 continue$c6 r0=r0*1.0e-3/nnode$c6 do 1002 n=1,nnode$c6 if (coorr(1,n).lt.r0)

28、coorr(1,n)=r0$c0 1002 continue$c6 epsilon=prmt(1)$c6 rho=prmt(2) 空一行stifdist=+u/x;u/x*ep*r+u/y;u/y*ep/r+u/z;u/z*ep*r 對應弱形式左端體積分項 空一行l(wèi)oad=+u*rou*r 對應弱形式右端體積分項 空一行end1.3 瞬態(tài)問題的“弱”形式 對于瞬態(tài)問題，求解其弱形式的過程與求解穩(wěn)態(tài)問題弱形式的過程基本相同，雖然瞬態(tài)問題的未知函數(shù)隨時間變化，在化為積分形式時，我們只是對空間變量積分。在pFEPG系統(tǒng)中填寫PDE(VDE)文件時，對于波動問題（即偏微分方程中有未知量關于時間的二階導

29、數(shù)），未知量關于時間的二階導數(shù)項用質量矩陣表示，未知量關于時間的一階導數(shù)項用阻尼矩陣表示；對于拋物問題（即偏微分方程中有未知量關于時間的一階導數(shù)），未知量關于時間的一階導數(shù)項用質量矩陣表示。 三維線彈性小變形動態(tài)問題的“弱”形式對于三維直角坐標問題，線彈性小變形波動問題基本方程如下：彈性體V域內任一點沿x,y,z方向的平衡方程為其中為表示體內任一點的應力狀態(tài)的6個應力分量，為正應力，為剪應力，為單位體積的體積力在x,y,z方向的分量，為密度常數(shù)，為密度常數(shù)乘以阻尼系數(shù)。在微小位移和微小變形的情況下，略去位移導數(shù)的高次冪，則應變和位移間的幾何方程和本構方程與本章節(jié)中的相同。彈性體的全部邊界為。一

30、部分邊界上已知外力。稱為力的邊界條件，這部分邊界用表示；另一部分邊界上彈性體的位移已知，稱為幾何邊界條件或位移邊界條件，這部分邊界用表示。這兩部分邊界構成彈性體的全部邊界，即力的邊界條件為：其中，為彈性體在邊界上單位面積的內力。設邊界外法線為N，其方向余弦為，則邊界上彈性體的內力可由下式確定：位移邊界條件為：其中，為彈性體在邊界上已知的位移。運用虛位移原理求位移由平衡方程，得其中，為三個方向的虛位移。由本章節(jié)中的分部積分公式將其化為弱形式，利用幾何方程可得上式表示為向量形式如下其中，為應力張量，為應變張量，為位移向量，為體積力向量，為面力向量（詳細的請見本章節(jié)）將本構方程代入上式，可得 將力和

31、位移邊界條件代入上式（注意，對于位移邊界條件，虛位移為0），可得由上面的弱形式，填寫的pFEPG系統(tǒng)文件如下：VDE文件（右邊文字是注釋）disp u v wcoor x y zfunc exx eyy ezz eyz exz exymate pe pv fu fv fw rou alpha pe,pv,rou,alpha分別對應弱形式中的, fu,fv,fw分別對應弱形式中的shap %1 %2gaus %3mass %1 rou 質量矩陣表示弱形式中未知量關于時間的二階導數(shù)項的體積分項damp %1 rou*alpha 阻尼矩陣表示弱形式中未知量關于時間的一階導數(shù)項的體積分項vect u

32、u v wvect f fu fv fwvect ev exx eyy ezzvect eg eyz exz exyvect dg (0.5-pv) (0.5-pv) (0.5-pv)matr dv 3 3(1.-pv) pv pvpv (1.-pv) pvpv pv (1.-pv) 空一行$c6 fact = pe/(1.+pv)/(1.-2.*pv) 空一行funcexx=+u/x 空一行eyy=+v/y 空一行ezz=+w/z 空一行eyz=+v/z+w/y 空一行exz=+w/x+u/z 空一行exy=+u/y+v/x 空一行stifdist=+ev_i;ev_j*dv_i_j*fac

33、t+eg_i;eg_i*dg_i*fact 對應弱形式左端體積分項,不包括關于時間導數(shù)的體積分項 空一行l(wèi)oad=+u_i*f_i 對應弱形式右端體積分項 空一行endFBC文件（右邊文字是注釋）disp u,v,w,coor x,y,mate fx fy fz 0.0;0.0;1.0d2; fx,fy,fz分別對應弱形式中的shap %1 %2gaus %3 空一行stifdist=+u;u*0.0 對應弱形式左端邊界積分項，例子中沒有，所以系數(shù)取0.0 空一行l(wèi)oad=+u*fx+v*fy+w*fz 對應弱形式右端邊界積分項 空一行end 三維瞬態(tài)熱傳導問題的“弱”形式對于三維直角坐標問題

34、，對于物體域瞬態(tài)熱傳導問題方程如下：其中，為溫度，是物體的沿x,y,z方向的導熱系數(shù)，是物體的密度，c是物體的比熱，是物體內熱源密度。物體的全部邊界為。一部分邊界上溫度已知這里邊界上溫度的給定值 另一部分邊界上熱流已知其中，是邊界上熱流的給定值，分別為邊界表面外法線在方向的方向余弦。在本例中，這兩部分邊界構成了物體的全部邊界，即運用虛位移原理求溫度利用瞬態(tài)熱傳導方程可得其積分形式為：其中，為溫度的虛位移。由本章節(jié)中的分部積分公式將其化為弱形式，可得將邊界條件代入上式（注意，對于已知溫度邊界條件，虛位移為0），可得由上面的弱形式，填寫的pFEPG系統(tǒng)文件如下：PDE文件（右邊文字是注釋）disp

35、 u,coor x,y,zshap %1 %2gaus %3mate rou ec ekx eky ekz eq rou,ec,eq分別對應弱形式中的ekx,eky,ekz分別對應弱形式中的mass %1 rou*ec 質量矩陣表示弱形式中未知量關于時間的一階導數(shù)項的體積分項 空一行stifdist=+u/x;u/x*ekx+u/y;u/y*eky+u/z;u/z*ekz 對應弱形式左端體積分項,不包括關于時間導數(shù)的體積分項 空一行l(wèi)oad=+u*rou*eq 對應弱形式右端體積分項 空一行endFBC文件（右邊文字是注釋）disp ucoor x,yshap %1 %2gaus %3mate

36、 eb 3.0e1; eb對應弱形式中的 空一行stifdist=+u;u*0.0 對應弱形式左端邊界積分項，例子中沒有，所以系數(shù)取0.0 空一行l(wèi)oad=+u*eb 對應弱形式右端邊界積分項 空一行end 二維不可壓縮流體瞬態(tài)Navier_Stokes方程的“弱”形式對于二維直角坐標問題，不可壓縮流體瞬態(tài)Navier_Stokes方程如下：對于流體域內，連續(xù)性方程為 其中，為的速度分量。動量方程為其中，為流體壓力，分別為流體域內方向上的單位外力；為流體密度常數(shù)，為粘性系數(shù)，為一常數(shù)。關于邊界條件，我們將在下面推導弱形式時提出，便易理解。利用虛位移原理求解首先，利用連續(xù)性方程可得：其中，為的虛

37、位移。利用動量方程可得：其中，分別為的虛位移。由本章節(jié)中的分部積分公式將其化為弱形式，可得 其中，為流體的全部邊界，將上面的三個弱形式合并可得：下面提出邊界條件，在流體進口處，邊界上的速度向量已知，即此時，虛位移為0。在固壁上，粘附邊界上速度向量為0，即此時，虛位移也為0；在流體出口處()， 將邊界條件代入上面的弱形式，可得由上面的弱形式，填寫的pFEPG系統(tǒng)文件如下：VDE文件（右邊文字是注釋）disp u,v,p coor x,ycoef un vn pn 導入上一迭代步的u,v,p,用于對非線性項的處理func div mate rou emu fx fy rou,emu,fx,fy分別

38、對應弱形式中的shap %1 %2gaus %3mass %1 rou rou 0.0 質量矩陣表示弱形式中未知量關于時間的一階導數(shù)項的體積分項vect x x yvect f fx fyvect u u vvect un un vn 空一行funcdiv=+u_i/x_i 空一行stifdist= +u_i/x_j;u_i/x_j*emu -div;p-p;div+u/x_j*un_j*rou+v/x_j*un_j*rou 對應弱形式左端體積分項,不包括關于時間導數(shù)的體積分項 空一行l(wèi)oad=+u_i*f_I*rou 對應弱形式右端體積分項 空一行end這個問題屬于非線性問題，需要線性化才能

39、計算，在這里我們用的是簡單迭代法來進行線性化的。1.4 求解解梯度的最小二乘法在實際問題中，常存在求得解后，需要算出解梯度物理量的情況，如固體力學中，求出位移，如何求應力；溫度場中，求出溫度，如何求熱流密度；電磁場中，求出電勢，如何求電場強度；等等。我們采用最小二乘法來求解的。 已知位移求應力對于三維線彈性小變形問題，存在已知位移求應力的問題，我們從線彈性本構方程出發(fā)，采用最小二乘法來求解。1. 線彈性本構方程 可將上式簡寫為 2. 用最小二乘法求應力對于彈性體，由本構方程得如下的弱形式由上面的弱形式，在填寫pFEPG系統(tǒng)的VDE文件時，弱形式的右端項是由質量矩陣表示的，而不是剛度矩陣。文件如

40、下：VDE文件（右邊文字是注釋）disp sx sy sz syz sxz sxycoef u v wcoor x y zload fsx fsy fsz fsyz fsxz fsxy 對應弱形式右端體積分項mate pe pv fu fv fw pe,pv分別對應弱形式中的,fu,fv,fw分別對應弱形式中的shap %1 %2gaus %3mass %1 1.0 質量矩陣表示弱形式的左端體積分項vect fsv fsx fsy fszvect fsg fsyz fsxz fsxyvect fs sx sy sz syz sxz sxyvect fss fsx fsy fsz fsyz fsxz fsxyvect ev ex ey ezvect eg eyz exz exyvect dg (0.5-pv) (0.5-pv) (0.5-pv)matr dv 3 3(1.-pv) pv pvpv (1.-pv) pvpv pv (1.-pv) 空一行$c6 fact = pe/(1.+pv)/(1.-2.*pv) 空一行stif$cv ex=+u/x$cv ey=+v/y$cv ez=+w/z$cv eyz=+v/z+w/y$cv