第4章弹性力学有限元

1、4-1 引言引言4-1 引言引言一、弹性力学的研究对象：二维、三维结构。一、弹性力学的研究对象：二维、三维结构。如板壳结构、水坝等。如板壳结构、水坝等。 解析法（精确法）解析法（精确法）二、二、 弹性力学的分析方法弹性力学的分析方法 数值法（近似法）数值法（近似法） 解析法只能解少数形状规则、边界条件简解析法只能解少数形状规则、边界条件简单、荷载简单的结构。单、荷载简单的结构。 数值近似（如差分法）数值近似（如差分法） 近似法近似法 物理近似（如反弯点法等）物理近似（如反弯点法等）:抓住问题的主要力学特征作某些近似假设。抓住问题的主要力学特征作某些近似假设。 将数值近似与物理近似结合，并使用计

2、将数值近似与物理近似结合，并使用计算机作为计算工具，便产生了有限单元法。算机作为计算工具，便产生了有限单元法。 三、有限单元法（有限元法）三、有限单元法（有限元法） 是一种运用计算机求解工程和科学问是一种运用计算机求解工程和科学问题的近似数值方法。题的近似数值方法。 有限元法的分析思路（步骤）：有限元法的分析思路（步骤）：1、单元划分、单元划分 将连续体划分为有限个将连续体划分为有限个有限大小有限大小的单元。的单元。单元之间通过结点相连。单元之间通过结点相连。 悬臂深梁悬臂深梁： j k j 单元与单元间单元与单元间 kk 通过结点相连通过结点相连： 每个单元的尺寸为有限大（不是无限小）每个单

3、元的尺寸为有限大（不是无限小） 称为有限单元，即有限元。称为有限单元，即有限元。 假设单元内部分布的位移场、应力场假设单元内部分布的位移场、应力场或应力或应力位移混合场。位移混合场。 以后仅讨论假设位移场的情况以后仅讨论假设位移场的情况 称有称有限元位移法，简称有限元法。限元位移法，简称有限元法。 2、单元分析、单元分析 将单元内各种量（如任一点的位移、应力、将单元内各种量（如任一点的位移、应力、应变等）均表示为单元结点位移的函数。应变等）均表示为单元结点位移的函数。 然后求单元刚度矩阵。然后求单元刚度矩阵。3、将各单元集合成整体、将各单元集合成整体 即形成结点平衡方程即形成结点平衡方程K D

4、 D=F。K 、F的形成方式同矩阵位移法。的形成方式同矩阵位移法。4、解方程组、解方程组K D D=F，求，求D D。5、计算单元内各点的应力、应变等。并对、计算单元内各点的应力、应变等。并对计算结果进行整理分析。计算结果进行整理分析。 可见，有限元法的步骤与矩阵位移法可见，有限元法的步骤与矩阵位移法基本相同。基本相同。 不同处：矩阵位移法得精确解，而有不同处：矩阵位移法得精确解，而有限元法得限元法得近似解近似解。 将连续体离散为有限个单元的集合。将连续体离散为有限个单元的集合。一、单元形状一、单元形状 二维单元（用于平面问题）二维单元（用于平面问题） 三维单元（用于空间问题）三维单元（用于空

5、间问题） 平面问题平面问题等参数单元等参数单元（模拟曲线边界）（模拟曲线边界） 将实际结构划分为若干单元的分割线称将实际结构划分为若干单元的分割线称为为网格网格。一般，网格越密（即单元数越多），。一般，网格越密（即单元数越多），计算结果越精确，但计算量增大。计算结果越精确，但计算量增大。计算结果计算结果 精确值精确值 单元数单元数 可见，当网格加密到一定程度，对计可见，当网格加密到一定程度，对计算结果的精度提高已有限。算结果的精度提高已有限。 故不应盲目加密网格。故不应盲目加密网格。 对较规则的问题，可均匀划分网格。对较规则的问题，可均匀划分网格。 应力较大（如应力集中）处，应加密应力较大（如

6、应力集中）处，应加密网格。网格。1、任意一个单元的结点，必须同时也是相邻单、任意一个单元的结点，必须同时也是相邻单元的结点，而不能是相邻单元边界上的内点。元的结点，而不能是相邻单元边界上的内点。 下图中，应在虚线处增加一根杆件。下图中，应在虚线处增加一根杆件。 会造成总刚会造成总刚K中主元不占优。中主元不占优。3、可将集中力的作用点或分布荷载的突变、可将集中力的作用点或分布荷载的突变点设置为结点。点设置为结点。4、在结构的厚度或材料性质有突变处，应、在结构的厚度或材料性质有突变处，应把突变线作为单元的分界线。把突变线作为单元的分界线。 结点编号时，应使每个结点与相邻结点编号时，应使每个结点与相

7、邻结点的编号之差尽可能小，以便使结点的编号之差尽可能小，以便使K的的带宽尽可能小。带宽尽可能小。 五、结点与单元的自由度五、结点与单元的自由度 结点自由度结点自由度：平面（：平面（2个）、空间（个）、空间（3个）个） 单元自由度单元自由度=结点自由度结点自由度单元结点数单元结点数l假定一种连续函数来描述单元内任一点的位假定一种连续函数来描述单元内任一点的位移移 称单元位移函数（模式）。称单元位移函数（模式）。l位移函数选取的好坏直接影响计算的收敛性位移函数选取的好坏直接影响计算的收敛性与精度。与精度。l通常选取多项式作为位移函数。通常选取多项式作为位移函数。1、广义坐标法、广义坐标法 利用帕斯

8、卡利用帕斯卡(Pascal)三角形。三角形。 1 常数项常数项 x y 线性项线性项 x2 xy y2 二次项二次项 x3 x2y xy2 y3 三次项三次项 x4 x3y x2y2 xy3 y4 四次项四次项 （1）对于一维单元）对于一维单元21123( )nnu xxxx（2）对于二维单元）对于二维单元22123456221234562( , )( , )niniiiiiiiu x yxyxxyyyv x yxyxxyyy 应依次选取常数项、线性项、二次应依次选取常数项、线性项、二次项、项、。 并应注意坐标轴方向的对称性。并应注意坐标轴方向的对称性。【例【例4-1】选取图示二结点轴力元的位

9、移函数。】选取图示二结点轴力元的位移函数。ujuiu(x)xE,A,I,l解：解：因有因有2个结点位移，可取个结点位移，可取 u(x)= 1+ 2x 121,()ijiuuul将将 1、 2代入，并整理得：代入，并整理得：(0) , ( )ijuuu lu由由 得：得： 其中，其中， Ni 、 Nj 称为称为形函数形函数， N称为称为形函数矩阵形函数矩阵。( )1ijxxu xuull ( )( )+ eiiijjijjuu xN uN uNNNu形函数的性质：本结点上取值形函数的性质：本结点上取值1、其它结点上取值、其它结点上取值0。写成矩阵形式，为写成矩阵形式，为1, ijxxNNll u

10、i i j uj 1 Ni Nj 1 形函数形函数Ni表示表示:当当第第i个个结点位移分量发生结点位移分量发生单单位位移位位移,而其它结点位移分量为零时整个单元的而其它结点位移分量为零时整个单元的变形形状变形形状,故称为故称为“形状函数形状函数”,简称形函数简称形函数。1ixNl jxNl2、插值函数法、插值函数法一维单元：一维单元： u(x)=N1u1+ N2u2+ Nnun=S S Niui = N (e) 形函数矩阵：形函数矩阵： N = N1 N2 Nn (e) =u1 u2 un T 二维单元：二维单元： u(x,y)=N1u1+ N2u2+ Nnun=S S Niui v(x,y)

11、=N1v1+ N2v2+ Nnvn=S S Nivi写成矩阵形式为写成矩阵形式为 u= N (e)其中其中 ( , )( , )u x yuv x y 1212 0 0 00 0 0 nnNNNNNNN ( )T1122 ennuvuvuv二维单元：二维单元： u(x,y)=N1u1+ N2u2+ Nnun=S S Niui v(x,y)=N1v1+ N2v2+ Nnvn=S S Nivi写成矩阵形式为写成矩阵形式为 u= N (e)二、有限元解答的收敛性二、有限元解答的收敛性 为了保证解答的收敛性，须使假设的为了保证解答的收敛性，须使假设的位移函数满足位移函数满足：（1）包含常数项）包含常数

12、项 反映单元的刚体位移。反映单元的刚体位移。（2）包含线性项）包含线性项 反映单元常应变。反映单元常应变。（3）除在单元内连续外，还应在相邻单元）除在单元内连续外，还应在相邻单元的公共边界上连续。的公共边界上连续。 分分 C0连续（位移连续）连续（位移连续） 等。等。 C1连续（位移的一阶导数连续）连续（位移的一阶导数连续） 一般要求一般要求C0连续。连续。 上述条件（上述条件（1）、（）、（2）是收敛的必要）是收敛的必要条件，加上条件（条件，加上条件（3）就是收敛的充分条件。）就是收敛的充分条件。 把满足条件（把满足条件（1）、（）、（2）的单元称为）的单元称为完备单元，满足条件（完备单元，

13、满足条件（3）的单元称为协调）的单元称为协调单元。单元。 对于完备且协调的单元，具有单调收对于完备且协调的单元，具有单调收敛性。敛性。 y Fyi i Fxi s sy t txy Fyj s sx Fym Fxm Fxj m j x 对三结点三角形单元：对三结点三角形单元： (e) = ui vi uj vj um vm T F(e)=Fxi Fyi Fxj Fyj Fxm Fym T 单元分析：单元分析：A、确定单元内任一点的位移确定单元内任一点的位移u、应力应力e e、应变、应变s s，并建立它们与，并建立它们与 (e)间的关系。间的关系。 B、求求 k (e) ，即，即F(e)与与 (

14、e)之间的变之间的变换矩阵。换矩阵。对于平面问题：对于平面问题： u=u vT e e= e ex e ey g gxy T s s= s sx s sy t txy T对于空间问题：对于空间问题： u=u v wT e e= e ex e ey e ez g gxy g gyz g gzx T s s= s sx s sy s sz t txy t tyz t tzxT u= N (e)如对于三结点三角形单元：如对于三结点三角形单元： u=Niui+ Njuj+ Nmum v=Nivi+ Njvj+ Nmvm N = Ni 0 Nj 0 Nm 0 0 Ni 0 Nj 0 Nmxvyuyvxu

15、xyyxgee,e e=L u, L 微分算子矩阵微分算子矩阵。对于平面问题：对于平面问题：写成矩阵形式：写成矩阵形式： 0 e e x u e e= e e y = 0 g gxy vxyyx 0 e e x u e e= e e y = 0 g gxy vxyyx即即 e e = L u对于空间问题：对于空间问题：xuxeyvyezwzexvyuxygywzvyzgxwzuzxge e = L u:xyx 0 0 0 0 0 0 0 0 0yzyxzzL=u= uvw 将将u= N (e)代入代入 e e=L u，得：，得： e e=L N (e) =B (e) 其中其中 B= L N 称

16、为应变矩阵或几何矩阵。称为应变矩阵或几何矩阵。 e e= B (e)以平面应力为例：（以平面应力为例：（ m m为泊松比）为泊松比） E s sx= ( e ex+m em ey ) 1 m m2 2 E s sy= ( e ey+m em ex ) 1 m m2 2 E 1 1 m m t txy= g gxy 1 m m2 2 2 2写成矩阵形式写成矩阵形式 s s=D e e，则，则 1 m m 0 0 E E m m 1 0 1 0 D= 1m m2 2 1 1 m m 0 0 2 若为平面应变问题，将若为平面应变问题，将E 、m m 分别替换为：分别替换为： E E 1 1 m m

17、1m m2 2 2 2 D 弹性矩阵。弹性矩阵。 将将e e= B (e)代入代入 s s=D e e得：得： s s=S (e)其中其中 S = D B= D L N 称为应力转换矩阵。称为应力转换矩阵。 利用虚功原理推导利用虚功原理推导 设单元结点力与结点位移列阵为：设单元结点力与结点位移列阵为： F(e)=F1 F2 F3 Fn T (e)= 1 2 3 nT 给一虚位移给一虚位移 * *(e)，则单元内任一点的，则单元内任一点的 虚位移虚位移 u *= N * *(e) ， 虚应变虚应变 e e* *= B * *(e) 。 由虚功原理：由虚功原理： （外力虚功）（外力虚功） W=U

18、（虚变形能）（虚变形能） W= 1* F1+ 2* F2+ n* Fn =( * *(e) )T F(e) U=V(e ex* * s sx+ e ey* * s sy+ g gzx*t tzx)dxdydz = Ve e* *Ts s dxdydz = ( * *(e) )T V B T s s dxdydz 由由W=U，得：，得： F(e) = V B T s s dxdydz 将将 s s=S (e) = D B (e) 代入，代入， 得得 F(e) = k (e) (e) 其中，其中， k (e) = V B T D B dxdydz 单元上非结点荷载单元上非结点荷载 （静力等效）（静

19、力等效） 单元等效结点荷载单元等效结点荷载 FE(e)静力等效静力等效 ： 原荷载与原荷载与FE(e)在任何虚在任何虚位移上作的虚功都相等。位移上作的虚功都相等。 设有设有F=Fx Fy Fz T作用在单元内作用在单元内(含含单元边界上单元边界上)某点某点Cx，y，z。 设有虚位移设有虚位移u *= N * *(e) ，则，则C点虚点虚位移为位移为uC*= NC * *(e) 。其中。其中 NC 为为 N 在在C点的取值。点的取值。 FE(e)作虚功：作虚功： ( * *(e) )T FE(e) ； F作虚功：作虚功： uC*T F= ( * *(e) )T NC T F FE(e)= NC

20、T F 设有分布面荷载设有分布面荷载f=fx fy fz T作用在单元某一边界面积作用在单元某一边界面积A上。上。 取一微元取一微元dA上的分布力作为集中力，上的分布力作为集中力，即即 dF=fdA ，代入上式得：，代入上式得：FE(e)=A N T dF FE(e) = A N T fdA 设单元内有一分布体力设单元内有一分布体力 f=fx fy fz T 将一微元将一微元dV上的体力上的体力 dF= fdV作为一集中力，则有：作为一集中力，则有： FE(e)=V N T dF FE(e) = V N T fdV 一、结构离散化一、结构离散化二、单元分析：二、单元分析： 1、确定位移模式、确

21、定位移模式 u= N (e) 2、分析应变和应力、分析应变和应力 e e= B (e) s s=D B (e) 3、计算单刚、计算单刚 k (e) = V B T D B dxdydz 4、计算单元等效结点荷载、计算单元等效结点荷载 FE(e) = A N T fdA(分布面力分布面力)FE(e) = V N T fdV(分布体力分布体力)三、整体分析三、整体分析 结点平衡方程组结点平衡方程组 K D D=F=FJ+FE k (e) FE(e) 由直接刚度法集成由直接刚度法集成 “对号入座，同号累加对号入座，同号累加”四、引入支承条件四、引入支承条件 主主1副副0法或乘大数法。法或乘大数法。五

22、、解方程组五、解方程组 K D D=F，求，求D D六、计算单元应力六、计算单元应力 s s=D B (e) (e)由由D D中取出。中取出。 对计算结果进行分析整理。对计算结果进行分析整理。 q qi(Mi) q qj(Mj) E, I, l x vi(FQi) i v(x) j vj(FQj) y x 图示平面梁单元的结点位移和结点力为：图示平面梁单元的结点位移和结点力为： (e)=vi q qi vj q qj T F(e)=FQi Mi FQj MjT1、单元位移、单元位移 v(x)= 1 1 2 2x+ 3 3x2 + 4 4x3 1 1、 2 2、 3 3 、 4 4由边界条件：由边界条件： x=0，v=vi； x=0，q q=dv/dx=q qi； x= l ， v=vj； x= l ， q q=dv/dx=q qj。求得。求得。 1vi 2qi 3 3( vi vj)/l2 (2qiqj)/l 42(vi vj) /l2 (qiqj)/l2 将1、 2、 3 、4代入 v(x)= 1 2x+ 3x2 + 4x3中，并整理得： v(x)=Nivvi+Niqqi+ Njvvj+ Njqqj=N (e)223323222332321