【正文】 ??4 ??1 2 3 1 2 P 32 有限元程序設(shè)計的一般步驟 1. 算法描述和列式推導(dǎo)； 2. 框圖設(shè)計； 3. 代碼編寫； 4. 上機調(diào)試、考 核； 5. 編寫應(yīng)用說明； 6. 修改、補充、完善。 程序設(shè)計的一般要求：具備較齊全的功能、較強的通用性和可移植性、較好的可擴充性、良好的可讀性、足夠的可靠性、良好的自適應(yīng)性。 輸入數(shù)據(jù)及分類 1. 控制數(shù)據(jù)：結(jié)點總數(shù)、單元總數(shù)、約束總數(shù) 、荷載總數(shù)、 問題類型數(shù)等 2. 幾何數(shù)據(jù)：結(jié)點坐標(biāo)、單元信息（各單元的結(jié)點編號）、約束條件、單元類型數(shù)（ 彈性模量 E、 泊松比 ?、厚度 t 不同為一類） 3. 物性數(shù)據(jù)： 彈性模量 E、泊松比 ?、厚度 t 4. 荷載數(shù)據(jù)：荷載類型（集中、分布）、位置、方向、大小等 三結(jié)點 三角形單元分析平面問題的主 體 程序 1. 程序框圖 （參見圖 ） 2. 結(jié)構(gòu)化程序設(shè)計方法 模塊化 —— 由 1 個主程序和若干子程序組成。子程序由通用性，采用可調(diào)數(shù)組；主程序采用動態(tài)數(shù)組存儲技術(shù) 3. 動態(tài)數(shù)組存儲技術(shù) （ 1）按整型和實型定義兩個大型共享數(shù)組，如 A(1000000)、 M(1000000)； （ 2）設(shè)計動態(tài)數(shù)組表：將各子程序中的變界（可調(diào)）數(shù)組按各自實際需要的大小分配一維數(shù)組的空間； （ 3）動態(tài)數(shù)組覆蓋技術(shù)：全部或部分覆蓋。 33 圖 平面問題主體分析程序框圖 讀入控制數(shù)據(jù) 開始 讀入幾何 、物性、荷載數(shù)據(jù) 平面應(yīng)力 /應(yīng)變 ? E=E0/(1?02), ?=?0/(1?0) E=E0, ?=?0 計算單元剛度元素 疊加到整體剛度矩陣中 計算等效結(jié)點荷載 引入位移約束條件 解線性代數(shù)方程得結(jié)點位移 計算單元應(yīng)力、反力等 輸出結(jié)果 結(jié)束 34 平面桿系結(jié)構(gòu)的有限元 等截面直梁單元（忽略剪切變形） (1) 基本方程 圖 受任意荷載的等截面直梁 幾何關(guān)系： 22dd,dd xvxux ??? ?? 寫成矩陣系數(shù)， }]{[}{ fLx ???????? ??? ， {f}=[ u v]T， ??????????????22dd00dd][xxL 內(nèi)力 位移關(guān)系： 22ddddxvEIEIMxuEAEANx??????? }]{[}{ ?? DMN ???????? ??????? EIEAD 0 0][ 平衡關(guān)系： x y p(x) q(x) 35 )(dddddddd)(dd223322xqxMxQxvEIxMQxpxuEA???????? 邊界條件： ? ?? ? 00:。,)(00:。,)(00:。,)(??????????????????MQMMMvMMvvvdxdvvvssssss，自由例如，簡支例如，固支例如 ?? (2) 離散 將一平面桿件結(jié)構(gòu)離散為 ne個單元， n 個結(jié)點，基本未知量為結(jié)點位移 {?}=[u1 v1 ?1 u2 v2 ?2 … un vn ?n]T 對應(yīng)的結(jié)點力向量為 {F}=[X1 Y1 M1 X2 Y2 M2 … Xn Yn Mn]T (3) 位移模式 局部坐標(biāo)下的單元結(jié)點位移向量和單元結(jié)點力向量： 圖 等截面直梁單元 TjjjiiieTjjjiiieMQNMQNFvuvu][}{][}{?? ??? 由梁的平衡方程可知，在結(jié)點荷載作用下梁的軸向位移沿梁軸呈線性分布 （因只有桿端作用軸向荷載） ，橫向位移呈三次曲線分布 （因彎矩沿桿軸成線性分布） ，故假設(shè) u=?1+?2 x v=?3+?4 x+?5 x2+?6 x3將結(jié)點 i, j 的位移代入可求出 6 個待定系數(shù) ?1~ ?6。 x y x’ y’ vi ui vj uj ?i ?j ? 36 將單元位移寫成結(jié)點位移的顯式，有 u=N1ui+N4 uj v= N2vi+ N3?i+ N5vj+ N6?j 23263322542323332221,23,2,231,1lxlxNlxlxNlxNlxlxxNlxlxNlxN?????????????? 位移模式的矩陣表示 eNvuf }]{[}{ ????????? ??????? 6532 41 00 0000][ NNNN NNN u、 v 獨立插值，但 ?不獨立插值，故要求 C1連續(xù)：不僅 u、 v 本身連續(xù)， v 的一階導(dǎo)數(shù)也要連續(xù)。 (4) 單元應(yīng)變矩陣和應(yīng)力矩陣 將位移模式 {f}=[N]{?}e代入幾何方程，得單元應(yīng)變?yōu)椋? ex BfL }]{[}]{[}{ ???? ????????? ???????????????22dd00d d][],][[][xxLNLB 將 {?}=[B]{?}e代入物理方程，可得 {?}=[D][B]{?}e, 其中 ??????? EIEAD 0 0][ (5) 單元剛度矩陣 局部坐標(biāo)下的單元剛度矩陣： ?? l T xBDBk d]][[][][ 37 單元剛度方程 為 ： {F}e =[k]{?}e [k]中的每一元素 krs稱為單元剛度系數(shù)，其物理意義是：第 s 個桿端位移為單位位移，而其他桿端位移為零時所引起的第 r 個桿端力。單元剛度矩陣是一個對稱矩陣，即 krs=ksr，這可由反力互等定理證明；單元剛度矩陣還是一個奇異矩陣，這是由于單元中包含剛體位移。 ???????????????????????????????????????lEIlEIlEIlEIlEIlEIlEIlEIlEAlEAlEIlEIlEIlEIlEIlEIlEIlEIlEAlEAk460260612061200000260460612061200000][222323222323 局部向整體坐標(biāo)的變換： {?’}e=[T]T{?}e {F’}e =[T]T{F}e [k’] =[T]T[k][T] 式中， [T]為一 6?6 的坐標(biāo)變換矩陣。假設(shè)自整體坐標(biāo) Ox’軸沿轉(zhuǎn)角正方向（即 順 時針 方向）轉(zhuǎn)到局部坐標(biāo) Ox 軸的角度為 ?，則 ???????????????????1000c o ss i n0s i nc o s][,0 0][ ????tttT {F’}=[k’]{?’}e 坐標(biāo)變換矩陣 [T]是一個正交矩陣，即 [T]?1=[T]T，故結(jié)點位移和結(jié)點力由整體向局部坐標(biāo)的變換式為 {?}e=[T]{?’}e， {F}e=[T]{F’}e 考慮剪切變形的梁單元 基本假定：垂直于梁軸線的截面在變形后仍保持為平面，但不再垂直于變形后的 38 軸線。 撓度由兩部分組成：彎曲變形部分和剪切變形部分，即 v=vb+vs； 桿軸線斜率： dv/dx=?+?（ ?：截面的轉(zhuǎn)角，由彎曲變形引起； ?：剪切變形 ，截面與桿軸夾角的改變 ） 圖 發(fā)生彎曲和剪切變形的直梁微段 (1) 在經(jīng)典梁單元基礎(chǔ)上引入剪切變形 的 影響 vb用三次式插值 ，參見上一節(jié)。 vs用線性插值： vs = N7vsi+ N8vsj lxNlxN ??? 87 ,1 。 由此可導(dǎo)得單元剛度矩陣為 ??????????????????????????????????????????????????)1()4()1(60)1()2()1(60)1(120)1(6)1(12000)1()4)1(60)1(120][2232323?????????????lEIlEIlEIlEIlEIlEIlEIlEAlEAlEIlEIlEIlEAk（對稱 式中，212GAlEIk??為剪切影響系數(shù)， k 為截面剪應(yīng)力不均勻分布修正系數(shù)，對矩形截面 k=，圓形截面 k=10/9。 (2) Timoshenko 梁單元 對 u、 v、 ?均獨立插值： u=N1ui+N2 u2， v= N1v1+ N2v2， ?=N1?1+ N2?2 式中 ?：截面的轉(zhuǎn)角， dv/dx：桿軸線的斜率。 2,)(2),1(21),1(21 2121 xxxl xxNN cc ???????? ??? dx dx ? ? 39 ???????????????????????????????????????????????????????????320620102100000320100100100000000010000][][][222lllllllklGAIAIAIAlEIkkksb對稱對稱 當(dāng) l/h??時只能得到零解，即將出現(xiàn)剪 切自鎖（ Shear locking）。這是由于 v、 ?同階插值 （實際上放大了剪切應(yīng)變項的量級） ，故 dv/dx 與 ?不同階，從而導(dǎo)致剪應(yīng)變?=dv/dx????不能處處滿足（ dv/dx 為常數(shù)， ?為一次式），除非 ?=常數(shù)，意味著梁不能發(fā)生彎曲。 解決方法： 減縮積分 —— 數(shù)值積分時采用比精確積分要求少的積分點數(shù)，例如對兩結(jié)點梁單元采用一點積分；減縮積分后 [ks]中的 l2/3 和 l2/6 項均改為 l2/4。 假設(shè)剪切應(yīng)變 —— 對剪應(yīng)變 ?另行假定插值形式； 替代插值函數(shù) —— 計算剪應(yīng)變時，對 ?采用低一階的插值函數(shù)，如兩結(jié)點梁單元其插值函數(shù)為常數(shù) 1/2，即 2211 ??? NN ?? 。 兩結(jié)點 Timoshenko 梁單元包含橫向剛體位移（ v=c）和剛體轉(zhuǎn)動（ ?=dv/dx=c），包含常剪切應(yīng)變（ ?=0， dv/dx=?=c），但不包含常彎曲應(yīng)變的位移狀態(tài)（ ?=cx，v=）， 不能分析純彎問題（將伴隨剪切應(yīng)變），因而 是一種較低級的 C0連續(xù)型梁單元，實際計算中多采用三結(jié)點或更多結(jié)點的梁單元。 板彎曲的有限元 薄板小 撓度 彎曲問題 (1) 基本假定 1）直法線假定： ?z=0， ?zx=0， ?yz=0； 2）不計 ?z引起的變形（物理方程與平面應(yīng) 力 問題相同）； 3）小撓度假定：變形后中面無面內(nèi)位移。 40 2. . 基本方程 幾何方程： {?}=[L]w， TyxwywxwL ????????????????? 22222 2][ 　　 ， 式中 {?}為曲（扭）率向量。 應(yīng)變 ： {?}=[?x ?y ?xy]T= z{?}。 物理方程： {M}=[Mx My Mxy]T=[D]{?} 式中 [D]為板彎曲的彈性關(guān)系矩陣， ????????????????21000101][ 0???DD , )1(12 230 ??? EtD為