计算结构力学课程讲义

第1章 绪 论1.1 课程内容(1) 研究内容本课程主要研究工程结构计算机分析（数值分析）的常用方法有限单元法、加权残数（余量）法和边界单元法的基本概念、基本原理及其应用。(2) 参考书籍课程的主要参考书籍如下：唐锦春，孙炳楠，郭鼎康，计算结构力学，浙江大学出版社，1989丁皓江, 谢贻权, 何福保，弹性和塑性力学中的有限单元法，机械工业出版社，1989王勖成，有限单元法，清华大学出版社，2003王勖成，邵敏，有限单元法基本原理与数值方法，第二版，清华大学出版社，1997徐次达，固体力学加权残数法，同济大学出版社，1987孙炳楠，项玉寅，张永元，工程中边界单元法及其应用，浙江大学出版社，1991Bath, K. J. Finite Element Procedures, Prentice-Hall, Inc., 1996.Zienkiewicz, O. C., The Finite Element Method, 5th Edition, McGraw Hill, 2001.Brebbia, C.A., The Boundary Element Method for Engineers, Pentech Press, London, 1978.Chandrupatla, T. R., Belegundu, A.D. Introduction to Finite Elements in Engineering, Prentice-Hall, Inc., 2002.1.2 结构分析方法概述一个工程技术问题总可由一组基本方程（通常是微分方程）加一组边界条件描述，即由下式给出：基本方程：L(u)-p=0， V（域内）边界条件：B(u)-g=0， S（边界）式中L、B为算子，p、g为已知函数。工程技术问题的常用分析方法有：(1) 解析方法只适用于少数简单问题，即形状规则且外部作用（如外荷载）简单的结构分析问题。(2) 数值方法数值方法可分为区域型方法和边界型方法。常用的区域型方法包括有限差分法、加权残数法、里兹(Ritz)法（变分法）和有限单元法等，其中有限差分法是直接对基本微分方程进行离散，再对离散后的代数方程进行求解；后几种方法则是先建立基本方程（一般是微分方程）的等效积分表达式，再进行离散求解。边界型方法中最典型的是边界单元法。它是先将基本微分方程变换为等效的边界积分方程，再在边界上对其进行离散求解。例如，图1.1给出了一个受复杂横向荷载（分布荷载、集中力、集中力偶等）作用的两端固定变截面梁。为求梁的挠度和内力，可列出梁的基本方程和边界条件如下：lq(x)yxh(x), EI(x)PmABbh(x)图1.1 变截面单跨梁受横向荷载作用基本方程：L(u)-p=0，V（域内）EI(x)y = -M(x), 0xl.边界条件：B(u)-g=0， S（边界）(y)x=0或x=l=0，(y )x=0或x=l=0以下分别就采用加权残数法、里兹法（位移变分法）和有限单元法的基本原理进行讨论。(1) 加权残数法为求近似解，设试探函数代入基本方程和边界条件，得残值：RL=L(u)-p（域内），RB=B(u)-g（边界）迫使残值在某种平均意义（加权积分）上等于零，则有由此可得到关于待定系数ai的代数方程组，解方程可求得待定系数及解答的近似表达式，其中的试函数可以选择多项式、三角函数、样条函数等。(2) 里兹法（位移变分法）里兹法的理论依据是最小势能原理。该原理可表述为：给定外力作用下，满足几何条件的各种可能位移中，真实的位移使总势能取极值，据此有d(U+UR)=0假设满足位移边界条件的位移函数为：将其代入方程得到关于待定系数Ai的代数方程，解方程可得Ai。里兹法需要在整个计算区域上假设近似函数，很难适应形状（边界）较复杂或解答较难预测的问题。(3) 有限单元法有限单元法的理论依据是最小势能原理或其他形式的变分原理。该方法与里兹法的主要区别是不在整体计算区域上假设近似函数，而是先将连续的求解区域离散为一个由有限个单元组成并按一定方式相互连接的单元集合体，再以各单元连接结点处的场量（如位移量）作为基本未知量，在各单元内假设近似函数（通过结点未知量插值得到），从而将一个无限自由度问题简化为有限自由度问题。u分段假设试函数x图1.2 一维试函数的分段假设例如图1.2中的曲线是某个一维问题的目标函数曲线，若采用里兹法对整个区段假设一个近似的试函数，显然比较困难。但如果现对整个区域进行分段（如图中短线为分段线），再对各个区段假设试函数，则要简单和准确得多，如可将各区段均假设为二次函数。哟次可见，有限单元法可视为一种分片（或分块、分段）形式的变分法。虽然有限单元法的理论依据和里兹法是一致的，但采用了分片（或分块、分段）假设试函数的处理方法以后，使得该方法的具体实施变得简便易行，具有了优越的可操作性和更为明确的物理意义，也使得该方法具有了其他方法（如里兹法）所不具备的优点：1) 概念简单、明确，易为工程人员接受；也可建立严格的数学分析和证明；2) 适用性十分广泛，适应于各类复杂边界和不同外部作用的问题；3) 求解过程程序化，易于编程和计算机实现。1.3 课时安排课程的总体课时安排如下：有限单元法部分包括概论、进展；平面三角形、矩形、等参元；杆元、板元等，共约20个课时；加权残数法部分包括基本原理、方法分类，以及伽辽金（Galerkin）方法、最小二乘法的应用，共需约46个课时；边界单元法主要包括基本原理（以二维势问题为例）；梁弯曲和板弯曲问题，共需约46个课时。思考题1.1 区域型分析法和边界型分析法在对问题的基本方程和边界条件的处理上有何不同和相同点？试分别举例说明。1.2 里兹法和有限单元法的理论依据、基本未知量的选取、试函数的假设等方面有何异同点？1.3 与里兹法相比，有限单元法在解决复杂问题上的适应性更为广泛，你认为主要的原因在于那些方面？第2章 有限单元法2.1 概述2.1.1 发展概况有限单元法的发展概况：1943年R. Courant尝试应用三角形区域上定义的分片连续函数和最小势能原理解决St. Venant扭转问题，是较早的有限元思想的体现：R. Courant, Variational Methods for the Solution of Problems of Equilibrium and Vibrations, Bulletin of the American Mathematical Society, 49: 1-23, 19431956年M.J. Turner，R.W. Clough等将刚架矩阵位移法推广到弹性力学平面问题，开始了有限元的第一个成功尝试和应用；用直接刚度法建立单元刚度特性：M.J. Turner, R.W. Clough, H.C. Martin and L.T. Topp, Stiffness and deflection analysis of complex structures, J. Aeronaut. Sci., 25: 805-823, 19561960年Clough第一次提出“有限单元法（FEM）”的名称，沿用至今。Zienkiewicz等编写第一本有限元方面专著：O.C. Zienkiewicz and Y.K. Cheung, The Finite Element Method in Continuum and Structural Mechanics, McGraw-Hill, New York, 19651963-1964年发现该方法是基于变分原理的里兹法的另一种形式，确立其理论基础。我国冯康在同一时期独立提出并证明了该方法：Melosh证明了位移法就是基于最小势能原理的Rayleigh-Ritz法冯康，基于变分原理的差分格式，应用数学和计算数学，1965，2(4): 238-2621960至今：实际工程应用：平面空间板壳；静力动力、波动稳定；弹性塑性粘弹性、复合材料；固体流体、传热等连续介质力学；计算分析优化设计、与CAD技术结合。E.L. Wilson：编写了第一个公开的有限元软件SAP；通用有限元软件：SAP、ADINA、NASTRAN、ANSYS、ABAQUS、MIDAS等从半个多世纪以来有限单元法的萌芽、理论依据的证明和充实及其逐步的广泛应用可以看到，它的发展和计算机软硬件的发展基本上是同步的。如果没有计算机的强大软硬件支撑，有限单元法只有其微不足道的一点理论上的意义，而没有更为重要的实际应用的意义。2.1.2 有限单元法概念(1) 离散化离散化的过程是将连续体划分为有限数目、有限大小的单元的集合体。单元与单元之间只在指定点（即结点）连接，其他位置则一般保持连续即可。单元可以具有不同的形状，即单元外形可以不同；单元与单元之间可以有不同的连接方式，即单元的结点数目、位置可以不同。连续体rg图2.1 连续体离散为单元集合体示例(2) 单元分析对典型单元假设位移模式（由各结点位移插值），再分析单元的力学特性，建立单元的结点力与结点位移之间的关系，即单元刚度方程：Fe=kDe）并将各类荷载变换为作用在结点上的等效结点荷载。(3) 整体分析将各单元刚度方程集成整体结构的整体刚度方程：F=KD根据结点的平衡条件，得最终的有限元方程：KD=R求解该方程可得到未知的结点位移。(4) 再次单元分析求出各单元的应变和应力。2.2 弹性力学平面问题的矩阵描述2.2.1两类平面问题弹性力学的平面问题可分为平面应力问题和平面应变问题两类。实际上，所有的弹性力学问题都是空间问题。所谓平面问题，并不是说这个问题所分析的对象本身（如形状、荷载分布）是平面的，而是指该问题的形状、外部作用以及问题的解答（即由此产生的效应，如位移、应力等）只在平面内有变化，而沿着平面外就保持不变了。因此可以肯定地说，所谓的平面问题就是一个特殊的空间问题。那么，是不是一个问题的形状和外部作用（即已知的位移和应力边界条件）只在平面内发生变化，而沿着平面外保持不变了，这个问题就是平面问题呢？不是的，还必须附加其他条件，这一结论才能成立。这个附加条件就是该问题沿平面外的尺寸与平面内尺寸相比要么非常小（如无限短），要么非常大（如无限长）。如果符合前者条件，则弹性体内只存在平面内的应力，而平面外的应力均为零，故这类问题称为平面应力问题；如果符合后者条件，则弹性体内只存在平面内的位移或平面内的应变，而平面外的位移及应变均等于零，故这类问题称为平面应变问题。zxt/2Oyt/2y体积力表面力图2.2 平面应力问题示例xOyz图2.3 平面应变问题示例2.2.2 基本量和基本方程的矩阵表示(1) 基本量应力：s=sx sy txyT，应变：e=ex ey gxyT，位移：f=u vT；体积力：p=px pyT，表面力：q=qx qyT。(2) 基本方程几何方程：物理方程：s=D e，式中D：弹性矩阵。对平面应力问题，平衡方程的弱形式能量原理1) 虚功原理：弹性体在外力作用下处于平衡状态的充分和必要条件是，外力在满足变形协调推进的虚位移上所做的虚功等于其内力在相应的虚应变上所做的虚功，即2) 最小势能原理：满足几何条件的各种可能位移中，真实位移使P取驻值，即dP=d (U+UR)=0(3) 平面问题的常用单元：三结点三角形单元 六结点三角形单元 四结点矩形单元四结点任意四边形单元 八结点曲边四边形单元图2.4 平面问题的常用单元单元的加密方法可分为p型加密和h型加密。前者是保持单元的大小和形状不变，而提高单元的插值函数的阶数；后者是指采用同一类单元，但加密单元的网格，即减小单元的尺寸，增加离散单元的数目。2.3 三结点三角形单元分析2.3.1 离散将一平面结构离散为ne个单元，设结点综述为n个，则结构的基本未知量取为这n个结点的位移，即D=u1 v1 u2 v2 un vnTyui (Ui)ixvi (Vi)uj (Uj) jvj (Vj)um (Um)mvm (Vm)图2.5 三结点三角形单元2.3.2 单元位移模式单元结点力向量：Fe=Ui Vi Uj Vj Um Vm T单元结点位移向量：De=ui vi uj vj um vmT(1) 位移模式的建立所谓位移函数或称位移模式，是指利用单元的结点位移将单元内任一点的位移用坐标的函数表示出来。对于三结点三角形单元，假设其位移模式是坐标的线性函数，则有u=a1+a2 x+a3 yv=a4+a5 x+a6y系数a1 a6由6个结点位移分量（ui、vi、uj、vj、um、vm）确定。将位移模式写成结点位移的显式：u= Niui+ Nj uj+ Nmumv= Nivi+ Nj vj+ NmvmNi、Nj、Nm称为形状函数（形函数）或插值函数，其中ai、bi、ci是分母行列式第1行各元素的代数余子式，展开后可表示为ai=xjym- xmyj，bi=yj-ym，ci= -xj+xm（i, j, m轮换） jyimxP j1imNi(x,y)图2.6 形函数的物理意义(2) 形函数性质(a) (Ni)i=1，(Ni)j=0，(Ni)m=0；(b) 单元内任一点：Ni+Nj+Nm=1。(3) 位移模式的矩阵表示其中N为形函数矩形，且。2.3.3 单元应变矩阵和应力矩阵将位移模式f=NDe代入几何方程，得单元应变为：e=BDeB称为单元应变矩阵，B= Bi Bj Bm，而各子块为。可见B为一常量矩阵，故三结点三角形单元为一常应变单元。将e=BDe代入物理方程，可得s=SDe=DBDeS称为单元应力矩阵，S= S1 S2 S3，对平面应力问题。可见S也是一常量矩阵，故三结点三角形单元为一常应力单元。2.3.4有限元方程的建立将连续体近似看作为由一系列只在结点相连的单元组装而成的集合体，并且对各个单元均建立了其位移模式，那么该集合体的最终位移法方程，即最终的有限元方程可由变分原理（此时为最小势能原理）建立。位移模式确定后，离散结构的可能位移就是由不同结点位移决定的分片函数。各可能位移函数中，真实的位移函数使离散结构的总势能取驻值（处于稳定平衡状态时为最小值），即dP=d (U+UR)=0令：De=GD这里G是一62n的位置矩阵，表示单元结点位移De在整体结点位移向量D中的位置；ke为单元刚度矩阵；Re单元等效结点荷载向量。于是有令：这里K为整体刚度矩阵；R为整体等效结点荷载向量，则有：由势能驻值原理dP=0可得或2.3.5 单元刚度矩阵(1) 单刚列式其中，(2) 单刚性质单元刚度元素kpq的物理意义为：第q个结点位移分量为单位位移（其它结点位移等于0），所引起的第p个结点力分量。例如要获得元素k26的直观解释，可先将单元的所有6个结点位移约束住，即在每个结点处沿两个方向同时施加链杆约束，如图2.7所示。ymxk26ij1图2.7 单元刚度元素k26的物理意义单元刚度矩阵k的性质： 1) 对称性： kpq = kqp； 2) 奇异性：因具有刚体位移； 3) 每行（列）元素之和为零；例如，第6列元素的意义：当第六个结点位移等于1（D6=1），其他结点位移等于0时，在各结点位移方向施加的结点力的大小，即图2.7中6个附加链杆约束上的约束力。因单元在这些结点力作用下处于平衡，故有：k16+ k36+ k56=0；k26+ k46+ k66=0。4) 主对角元素恒大于零：kqp0； 5) k取决于单元的形状、方位和弹性常数，与所在位置（即平移或np 转动）无关。(3) 推导单刚的另一种方法由单元的平衡条件导出物理方程几何方程位移模式虚功方程f =Ndee=Bdes=Sde ，S= DBFe=kd e，k= BT D BtA结点位移 位移 应变 应力 结点力 de f e s Fe图2.8 单元分析中各物理量的联系图假设单元发生了虚位移，则结点力在结点虚位移上做的虚功=单元应力在虚应变上的虚功（虚变形能）：令：又因dDe是任意的，故有Fe=keDe2.3.6 单元等效结点荷载在将连续体用近似的单元集合体代替以后，单元集合体中的每个结点可认为一方面受到外荷载的作用，另一方面受到连接该结点的各个单元的对其产生的约束力（即结点力）的作用。各结点在这两组力的作用下总是处于平衡的。但是，如果外荷载不是直接作用在结点上，而是作用在单元内部，那么怎么建立结点的平衡关系呢？一种有效方法是先建立单元的等效结点荷载，再对结点建立平衡关系。yYiixXiYj jXjYmmXmPxPy建立单元的等效结点荷载，或称将非结点荷载需等效移置到结点上，该等效移置可利用虚功原理进行。图2.9 单元的等效结点荷载分量2.3.7 整体刚度矩阵(1) 建立有限元方程的方法：1) 由最小势能原理建立：K=(Ge)Tk e Ge2) 由结点平衡条件建立(2) 整体刚度矩阵K的集成方法：对号入座。K=(Ge)Tk e Ge，其中(Ge)T决定kpqe在K中的行，Ge决定在K中的列位置。(3) 整体刚度矩阵的性质K元素的物理意义 Kpq：结构第q个结点位移分量为单位位移（Dq=1，其它结点位移=0），在第p个结点位移方向所施加的结点力的大小。 如 K25为结构第5个结点位移分量为单位位移，即D5=1，其它结点位移均为零时，在第2个结点位移方向所需施加的结点力的大小。K的性质： 1) 对称性： Kpq = Kqp；主对角元恒大于零； 2) 稀疏矩阵，且一般呈带状分布；例如，平面问题最大半带宽=2(单元结点号之差最大值+1)。 3) 引入位移边界条件前为奇异矩阵，引入后为正定矩阵。2.3.8 位移边界条件的引入以上集成的总体整体刚度矩阵并不包含位移边界条件的信息，这种刚度矩阵常称为原始刚度矩阵，它是一个奇异矩阵。要使得最终的有限元方程可解，必须引入唯一边界条件，排除结构的刚体位移。引入位移边界条件的常用方法有以下几种：(1) 直接引入法（矩阵缩小法）将KD=R改写为如下子块形式：其中Da、Db分别为未知结点位移向量和已知结点位移向量，Kaa、Kab、Kba、Kbb、Ra、Rb分别为与之对应的刚度矩阵和荷载向量子块。将上述方程的上半部分（对应已知位移）展开，得Kaa Da=Ra-Kab Db由该方程可解出未知结点位移。该方法将改变原始刚度矩阵的阶数及元素的顺序，不利于程序实现，因此在计算机编程中一般很少采用。(2) 对角元素改1法（零位移边界）设结构的总自由度数（即结点位移总数）为N，若第i个结点位移为零（即Di=0），则将K中对角元素Kii改为1，而第i行和第i列的其他元素改为0，荷载向量R中的第i个元素也改为0。其实质是将原有限方程中的第i个方程用方程Di=0代替，而其他方程中与Di对应的系数也改为0，表明Di对其他方程没有影响，同时保证了修改后的刚度矩阵仍具有对称性。(3) 乘大数法若第i个结点位移已知，即，则将整体刚度矩阵K中的对角元素Kii改为a Kii，其中a为一大数（如1020），而荷载向量R中的第i个元素Ri改为，原方程成为以下形式其实质是将原有限方程中的第i个方程用的以下近似方程代替：将上述方程各项同除以大数a，除第i项及方程右端项外，其他各项均趋于0，故等价于。2.3.9 位移模式与解答的收敛准则(1) 有限元解答的收敛准则为使有限元解答能够收敛于精确解，单元位移模式应满足以下条件：1) 位移模式必须包含单元的刚体位移；对弹性力学平面问题，其刚体位移表达式为：u=u0-wy，v=v0+wx因此，平面问题的位移模式必要包含上述刚体位移表达式中的各项，才能保证最终解答的收敛性。2) 位移模式必须包含单元的常量应变；条件(1)、(2)合并起来可称为完备性要求。对平面问题来说，就是要求位移模式必须包含常数项和完整的一次项。完备性条件是解答收敛的必要条件。3) 位移模式应尽可能保证位移的连续性。该条件实际上就是协调性条件，但一般情况下并不是一个必要性条件。如果位移模式同时满足上述完备性和协调性条件，那么就组成了解答收敛的充分条件。对一般的弹性力学平面或空间问题，只需要求单元内部以及相邻单元的公共边界上的位移本身连续，即位移模式具有C0连续性。对有些问题，可以放松对协调性的要求，只要通过分片试验，那么也能保证解答的收敛性。三结点三角形单元的位移模式既满足完备性，又满足协调性的要求（在单元边界上呈线性分布，可由两个结点位移唯一确定），是一种协调单元。证明如下：单元内：单值连续；yx(1)ijmp(2)相邻单元之间：uij(1)= uij(2)？ vij(1)= vij(2)？ij边的方程：y=ax+b，则uij=a1+ a2x+a3(ax+b)= cx+duij(1)、uij(2)均为坐标的线性函数，故可由i、j两点的结点位移唯一确定。图2.10 两相邻单元(2) 多项式位移模式的一般选择规则位移模式应与单元局部坐标的选取无关，即满足所谓的几何各向同性。对于平面问题，位移模式中的x和y的各阶项应保持对称，即有了xmyn项，则应同时具有xnym项。为保证位移模式关于x和y坐标的对称性，通常从以下的Pascal三角形中选取多项式位移模式的各项：1x yx2 xy y2 x3 x2y x y2 y3 x4 x3y x2y2 xy3 y4图2.11 Pascal三角形根据最小势能原理建立的有限元，是以结点位移作为基本未知量，这种有限单元称为位移元。由位移元得到的近似解答总体上是精确解的一个下限。2.4 高精度的三角形单元、矩形单元2.4.1 高精度三角形单元(1) 六结点三角形单元（二次单元T6）位移模式取坐标的完整二次式：u=a1+a2 x+a3 y+a4 x2+a5 xy+a6 y2v=b1+b 2 x+b 3 y+b 4 x2+b5 xy+b6 y2该位移模式包含了常数项和完整的一次项，满足完备性要求；在单元的边界上位移呈二次抛物线分布，可由三个结点位移唯一确定，故又满足协调性的要求，是一种协调单元。(2) 十结点三角形单元（三次单元T10）位移模式取坐标的完整三次式：u=a1+a2 x+a3 y+a4 x2+a5 xy+a6 y2+a7x3+a8x2y+a9xy2+a10y3v=b1+b 2 x+b 3 y+b 4 x2+b5 xy+b6 y2+b7x3+b8x2y+b9xy2+b10y3该位移模式包含了坐标的完整一次式（常数项和纯一次项），满足完备性要求；在单元的边界上位移呈三次曲线分布，可由4个结点位移唯一确定，故又满足协调性的要求，是一种协调单元。ijmijm3122.12 高精度三角形单元(3) 面积坐标表示的三角形单元形函数在推导三角形单元的列式时，直接利用整体坐标系（为直角坐标）下的位移模式将使得运算十分繁琐和复杂。如果采用三角形单元内的一种局部坐标面积坐标作为自然坐标，则可以使列式推导大为简化。定义单元内任一点P的无量纲面积坐标（Li, Lj, Lm）为i(1,0,0)Pj(0,1,0)m(0,0,1)AiAjAmLi= Ai/ A （i, j, m）各种单元的形函数：3结点三角形单元（线性单元T3）：Ni=Li （i, j, m）6结点三角形单元（二次单元T6）：Ni=(2Li-1) Li （i, j, m）N1=4Lj Lm, （1, 2, 3；i, j, m）2.13 三个结点单元示意图面积坐标的特点：1) 三角形三个角点的面积坐标：i(1,0,0)、j(0,1,0)、m(0,0,1)三条边用面积坐标表示的方程为：jm边 Li=0mi边 Lj=0ij边 Lm=02) 三个面积坐标不独立，其相互关系为Li+ Lj+ Lm=1面积坐标与直角坐标之间的关系：或2.4.2 四结点矩形单元（R4单元）(1) 采用双线性的位移模式：u=a1+a2 x+a3 y+a4xy1yx234hxov=b1+b 2 x+b 3 y+b 4xy可保证位移模式的完备性和协调性。若写成形函数形式，则为2.14 四结点矩形单元这里的Ni（i=1, 2, 3, 4）可以先求出8个待定系数在获得，也可以根据形函数的性质直接构造，例如对N1，可设该表达式满足在结点2、3、4取值均为零的性质；再令Ni(1)=1，可得待定系数这里设矩形单元的边长各为2a、2b。(2) 局部坐标下的形函数设矩形长和宽各为2a和2b，在矩形形心为原点建立局部坐标（自然坐标）系xoh，它与直角坐标的变换式：x4h123x= -1x=1h= -1h=1o2.15 局部坐标系下的矩形单元则四个角点的局部坐标分别为1(-1, -1)，2(1, -1)，3(1,1)，4(-1,1)。位移模式为：Ni（i=1, 2, 3, 4）可根据形函数的性质直接构造出来： 或(3) 单元应变矩阵和应力矩阵单元应变为：e=BDe其中B=B1 B2 B3 B4，。单元应力为：s=SDe=DBDe其中S= S1 S2 S3 S4，Si=D Bi2.5 C0连续型单元形函数的构造在有限单元法中，根据结点布置方式的不同可将矩形单元分为两类，一类称为Lagrange矩形单元，另一类称为Serendipity矩形单元。前者在单元纵横网格线的交点上均布置结点，而后者仅在单元的边界线上布置结点，如图2.17所示。(a) Lagrange矩形单元 (b) Serendipity矩形单元2.16 两类矩形单元2.5.1 Lagrange矩形单元(1) 一维Lagrange插值多项式过n个结点（坐标分别为x1, x2, , xn）的一维Lagrange插值多项式为xx1x2x3xnxi12.17 一维Lagrange插值多项式l1(x)例如n=2，则有：。若令x1=0，x2= l，则。用该多项式作为形函数，可满足形函数的性质。(2) Lagrange矩形单元的形函数在矩形的各个网格交点上均布置结点，如水平方向r+1个，竖向p+1个。于是第I列J行结点i的相应形函数：Ni= NIJ= lI(r)(x) lJ(p)(h )其中，lI(r)(x)、lJ(p)(h )均为一维Lagrange插值多项式，Ni在结点i为1，在其他结点处为零；单元边界上的结点数=形函数阶数+1，能够保证边界位移的唯一性和协调性。这类单元包含较多的内部结点，增加了单元的自由度，而实践证明这些自由度的增加通常并不能有效提高单元的精度。对nn结点的Lagrange单元，其Pascal三角形中包含的项数如图2.19所示。(I,J)(r, p)(0, 0)11lI (x)lJ (h)1xnynxnynxx2x3yy2y3图2.18 Lagrange矩形单元及插值模式 图2.19 Pascal三角形包含项数2.5.2. Serendipity矩形单元(1) 结点布置特点这类单元只在其边界上布置结点，但不同边界上可布置有不同数目的结点。(2) 形函数的构造方法如果只有四个角点上布置结点（R4单元），则如果增加结点5使之成为5结点矩形单元（或称R5单元），则满足(N5) j=d5j （j=1, 2, 3, 4, 5）。原不满足，需对其进行修正才能满足：(N1)5=0，(N2)5=0，且保证N1，N2在除5结点外的其他4个结点取值不变。x4h12351011/21N5图2.20 五结点矩形单元 图2.21 形函数N1的修正对于8结点Serendipity矩形单元（R8单元），利用上述方法，可得到前4个结点对应的形函数分别为或写成统一表达式：后4个形函数分别为对于p次单元的边界内结点，其相应的形函数为x（或h）的p次和h（或x）的一次Lagrange多项式的乘积；而角结点的形函数为四结点单元的相应函数与各内结点形函数乘以一分数的差。利用该方法可以很方便地构造一些过渡单元的形函数。2.6 平面等参单元直接用形状规则的单元对复杂形状的结构进行离散绘遇到边界难于处理的问题，因此应设法采用适当的方法对规则单元进行坐标变换，使之转化为斜边或曲边的单元。有限元中最常用的变换方法是等参变换，即坐标变换采用与单元位移插值一致的形函数。通过等参变换的单元称为等参单元。2.6.1 四结点四边形等参元对于4结点任意四边形单元，若采用双线性的位移模式，则位移在单元斜边界上呈二次抛物线分布，由两个结点位移不能唯一确定。为此在单元内设法建立一个局部坐标（自然坐标），使得单元各边界的局部坐标分别为一常量（或1），这样，如果位移模式采用局部坐标的双线性函数，则能够满足公共边界位移的唯一性。如果单元边界的局部坐标分别为1，则在局部坐标下，单元的形态就是图示的4结点正方形单元，该单元称为基本单元或母单元，其形函数为。位移模式为：1yx234xhx4h123x= -1x=1h= -1h=1o(a) 基本单元 (b) 实际单元2.22 四结点等参数单元如果同时利用上述形函数作为局部（自然）坐标(x, h)向整体（直角）坐标(x, y)的变换函数，则可以建立两种坐标之间的映射关系如下：根据形函数的性质：(Ni) j=dij，SNi=1，容易验证上述变换式在四个结点处给出了节点的整体坐标，而在单元的四边上，其中一个局部坐标为1，另一局部坐标按线性变化，因而正好给出了整体坐标下四边的直线方程（由两个结点整体坐标唯一确定）。例如对12边，因此上述变换式正确地反映了局部（自然）与整体（直角）坐标之间的映射关系。经变换后的实际单元称为子单元。子单元的位移模式仍为：利用等参变换可以构造8结点、12结点、20结点等更高次的四边形曲边等参单元（参见图2.23）。hxo15247386yx15247386图2.23 8结点曲边四边形等参单元2.6.2 局部与整体坐标的微分和积分变换式根据复合函数的求导法则，写成矩阵形式，式中J称为雅可比（Jacobi）矩阵，对于具有m个结点的平面等参单元，若反过来用局部坐标表示整体坐标，则可对上式作反变换，式中，Jij*是J各元素Jji的代数余子式；J是J的行列式，称为雅可比（Jacobi）行列式。面积微分的变换：dA=dxdy=Jdxdh。2.6.3 单元刚度矩阵和单元等效结点荷载向量对于一具有m个结点的平面等参单元，其应变向量可写为e=BDe其中，单元应变矩阵：B=B1 B2 Bm。单元应力向量：s=SDe=DBDe单元刚度矩阵：单元等效结点荷载列阵：式中Rpe、Rqe分别为体积力和表面力引起的等效结点荷载，且对其中一个局部坐标（x或h）为常数的边界，线积分ds为或上述积分式通常采用Newton-Cotes或高斯（Gauss）数值积分求得。可见，等参单元的刚度矩阵、等效结点荷载列阵等的计算都在规则的母单元中进行，因此各种形式的积分运算都可以采用标准的数值积分方法进行，使得不同工程问题的有限元分析能够采用统一的通用化程序实现。2.6.4 等参变换的应用条件等参单元的应用条件是母单元与实际单元之间存在一一对应关系。具体到局部与整体坐标的变换式，要求变换矩阵即雅可比J非奇异，雅可比（Jacobi）行列式J0。为确保坐标变换一一对应，在单元划分时应避免：(1) 使任意两个结点退化为一个结点而使dx=0或dh=0；(2) 还应避免因单元过于歪斜而导致dx与dh发生共线。xh123x= -1x=1h= -1h=1o4xh123,4x= -1x=1h= -1h=1o(a) 结点退化 (b) 单元边界接近共线2.24 畸形等参数单元2.6.5 例题分析例2.1 图示悬臂平面结构，长2m，高1m，试用图示单元进行分析。12346(1)(2)(3)5(4)xyP jmi(1)(3) jmi(2)(4)(a) 基本单元 (b) 实际单元2.25 例2.1图（三角形单元）1. 三结点三角形单元分析(1) 离散化：划分为4个单元，共6个结点。(2) 依据图示的单元局部编号规则，4个单元的刚度矩阵均相同，为为简便起见，设泊松比n=0，于是有(3) 根据局部与整体结点编号的对应关系集成整体刚度矩阵：(4) 引进位移边界条件2. 四结点四边角形等参单元分析12346(1)(2)5xyPx4h1232.26 例2.1图（四边形单元）计算步骤：(1) 离散化：划分为2个单元，共6个结点。单元上下短边0.75m，长边1.25m。(2) 求，J，|J|，J-1；(3) 求B在各积分点的数值Big；(4) 利用高斯积分计算并形成k；(5) 集成K、P；(6) 引进位移边界条件。12P例2.2 牛腿受竖向集中力作用，且结点1、2发生已知的水平位移，求应力。2.27 例2.2图2.7 有限元程序实现有限元程序前处理程序：生成网格及数据文件主体分析程序：核心计算分析后处理程序：进行结果处理，生成可直接使用或查看的结果文件、图形文件。2.7.1 有限元程序设计的一般步骤1. 算法描述和列式推导；2. 框图设计；3. 代码编写；4. 上机调试、考核；5. 编写应用说明；6. 修改、补充、完善。程序设计的一般要求：具备较齐全的功能、较强的通用性和可移植性、较好的可扩充性、良好的可读性、足够的可靠性、良好的自适应性。2.7.2 输入数据及分类1. 控制数据：结点总数、单元总数、约束总数、荷载总数、问题类型数等2. 几何数据：结点坐标、单元信息（各单元的结点编号）、约束条件、单元类型数（弹性模量E、泊松比n、厚度t不同为一类）3. 物性数据：弹性模量E、泊松比n、厚度t4. 荷载数据：荷载类型（集中、分布）、位置、方向、大小等2.7.3 三结点三角形单元分析平面问题的主体程序1. 程序框图（参见图2.28）2. 结构化程序设计方法模块化由1个主程序和若干子程序组成。子程序由通用性，采用可调数组；主程序采用动态数组存储技术3. 动态数组存储技术（1）按整型和实型定义两个大型共享数组，如A(1000000)、M(1000000)；（2）设计动态数组表：将各子程序中的变界（可调）数组按各自实际需要的大小分配一维数组的空间；（3）动态数组覆盖技术：全部或部分覆盖。读入控制数据开始读入几何、物性、荷载数据平面应力/应变?E=E0/(1-n02),n=n0/(1-n0)E=E0,n=n0计算单元刚度元素叠加到整体刚度矩阵中计算等效结点荷载引入位移约束条件解线性代数方程得结点位移计算单元应力、反力等输出结果结束图2.28 平面问题主体分析程序框图2.8 平面杆系结构的有限元2.8.1 等截面直梁单元（忽略剪切变形）(1) 基本方程xyp(x)q(x)图2.29 受任意荷载的等截面直梁几何关系：写成矩阵系数，f= u vT，内力-位移关系：平衡关系：边界条件：(2) 离散将一平面杆件结构离散为ne个单元，n个结点，基本未知量为结点位移D=u1 v1 q1 u2 v2 q2 un vn qnT对应的结点力向量为F=X1 Y1 M1 X2 Y2 M2 Xn Yn MnT(3) 位移模式局部坐标下的单元结点位移向量和单元结点力向量：xyxyviuivjujqiqja图2.30 等截面直梁单元由梁的平衡方程可知，在结点荷载作用下梁的轴向位移沿梁轴呈线性分布（因只有杆端作用轴向荷载），横向位移呈三次曲线分布（因弯矩沿杆轴成线性分布），故假设u=a1+a2 xv=a3+a4 x+a5 x2+a6 x3将结点i, j的位移代入可求出6个待定系数a1 a6。将单元位移写成结点位移的显式，有u=N1ui+N4 ujv= N2vi+ N3qi+ N5vj+ N6qj位移模式的矩阵表示u、v独立插值，但q不独立插值，故要求C1连续：不仅u、v本身连续，v的一阶导数也要连续。(4) 单元应变矩阵和应力矩阵将位移模式f=NDe代入几何方程，得单元应变为：将e=BDe代入物理方程，可得S=DBDe,其中(5) 单元刚度矩阵局部坐标下的单元刚度矩阵：单元刚度方程为：Fe =kDek中的每一元素krs称为单元刚度系数，其物理意义是：第s个杆端位移为单位位移，而其他杆端位移为零时所引起的第r个杆端力。单元刚度矩阵是一个对称矩阵，即krs=ksr，这可由反力互等定理证明；单元刚度矩阵还是一个奇异矩阵，这是由于单元中包含刚体位移。局部向整体坐标的变换：De=TTDeFe =TTFek =TTkT式中，T为一66的坐标变换矩阵。假设自整体坐标Ox轴沿转角正方向（即顺时针方向）转到局部坐标Ox轴的角度为a，则F=kDe坐标变换矩阵T是一个正交矩阵，即T-1=TT，故结点位移和结点力由整体向局部坐标的变换式为De=TDe，Fe=TFe2.8.2 考虑剪切变形的梁单元基本假定：垂直于梁轴线的截面在变形后仍保持为平面，但不再垂直于变形后的轴线。挠度由两部分组成：弯曲变形部分和剪切变形部分，即v=vb+vs；