计算结构力学课程讲义

PAGEPAGE72第1章绪论1.1课程内容(1)研究内容本课程主要研究工程结构计算机分析（数值分析）的常用方法——有限单元法、加权残数（余量）法和边界单元法的基本概念、基本原理及其应用。(2)参考书籍课程的主要参考书籍如下：唐锦春，孙炳楠，郭鼎康，计算结构力学，浙江大学出版社，1989丁皓江,谢贻权,何福保，弹性和塑性力学中的有限单元法，机械工业出版社，1989王勖成，有限单元法，清华大学出版社，2003王勖成，邵敏，有限单元法基本原理与数值方法，第二版，清华大学出版社，1997徐次达，固体力学加权残数法，同济大学出版社，1987孙炳楠，项玉寅，张永元，工程中边界单元法及其应用，浙江大学出版社，1991Bath,K.J.FiniteElementProcedures,Prentice-Hall,Inc.,1996.Zienkiewicz,O.C.,TheFiniteElementMethod,5thEdition,McGrawHill,2001.Brebbia,C.A.,TheBoundaryElementMethodforEngineers,PentechPress,London,1978.Chandrupatla,T.R.,Belegundu,A.D.IntroductiontoFiniteElementsinEngineering,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给出了一个受复杂横向荷载（分布荷载、集中力、集中力偶等）作用的两端固定变截面梁。为求梁的挠度和内力，可列出梁的基本方程和边界条件如下：llq(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（边界）迫使残值在某种平均意义（加权积分）上等于零，则有由此可得到关于待定系数i的代数方程组，解方程可求得待定系数及解答的近似表达式，其中的试函数可以选择多项式、三角函数、样条函数等。(2)里兹法（位移变分法）里兹法的理论依据是最小势能原理。该原理可表述为：给定外力作用下，满足几何条件的各种可能位移中，真实的位移使总势能取极值，据此有(U+UR假设满足位移边界条件的位移函数为：将其代入方程得到关于待定系数Ai的代数方程，解方程可得Ai。里兹法需要在整个计算区域上假设近似函数，很难适应形状（边界）较复杂或解答较难预测的问题。(3)有限单元法有限单元法的理论依据是最小势能原理或其他形式的变分原理。该方法与里兹法的主要区别是不在整体计算区域上假设近似函数，而是先将连续的求解区域离散为一个由有限个单元组成并按一定方式相互连接的单元集合体，再以各单元连接结点处的场量（如位移量）作为基本未知量，在各单元内假设近似函数（通过结点未知量插值得到），从而将一个无限自由度问题简化为有限自由度问题。uu分段假设试函数x图1.2一维试函数的分段假设例如图1.2中的曲线是某个一维问题的目标函数曲线，若采用里兹法对整个区段假设一个近似的试函数，显然比较困难。但如果现对整个区域进行分段（如图中短线为分段线），再对各个区段假设试函数，则要简单和准确得多，如可将各区段均假设为二次函数。哟次可见，有限单元法可视为一种分片（或分块、分段）形式的变分法。虽然有限单元法的理论依据和里兹法是一致的，但采用了分片（或分块、分段）假设试函数的处理方法以后，使得该方法的具体实施变得简便易行，具有了优越的可操作性和更为明确的物理意义，也使得该方法具有了其他方法（如里兹法）所不具备的优点：1)概念简单、明确，易为工程人员接受；也可建立严格的数学分析和证明；2)适用性十分广泛，适应于各类复杂边界和不同外部作用的问题；3)求解过程程序化，易于编程和计算机实现。1.3课时安排课程的总体课时安排如下：有限单元法部分包括概论、进展；平面三角形、矩形、等参元；杆元、板元等，共约20个课时；加权残数法部分包括基本原理、方法分类，以及伽辽金（Galerkin）方法、最小二乘法的应用，共需约4～6个课时；边界单元法主要包括基本原理（以二维势问题为例）；梁弯曲和板弯曲问题，共需约4～6个课时。思考题区域型分析法和边界型分析法在对问题的基本方程和边界条件的处理上有何不同和相同点？试分别举例说明。里兹法和有限单元法的理论依据、基本未知量的选取、试函数的假设等方面有何异同点？与里兹法相比，有限单元法在解决复杂问题上的适应性更为广泛，你认为主要的原因在于那些方面？

第2章有限单元法2.1概述2.1.1发展概况有限单元法的发展概况：1943年R.Courant尝试应用三角形区域上定义的分片连续函数和最小势能原理解决St.Venant扭转问题，是较早的有限元思想的体现：R.Courant,VariationalMethodsfortheSolutionofProblemsofEquilibriumandVibrations,BulletinoftheAmericanMathematicalSociety,49:1-23,19431956年M.J.Turner，R.W.Clough等将刚架矩阵位移法推广到弹性力学平面问题，开始了有限元的第一个成功尝试和应用；用直接刚度法建立单元刚度特性：M.J.Turner,R.W.Clough,H.C.MartinandL.T.Topp,Stiffnessanddeflectionanalysisofcomplexstructures,J.Aeronaut.Sci.,25:805-823,19561960年Clough第一次提出“有限单元法（FEM）”的名称，沿用至今。Zienkiewicz等——编写第一本有限元方面专著：O.C.ZienkiewiczandY.K.Cheung,TheFiniteElementMethodinContinuumandStructuralMechanics,McGraw-Hill,1963-1964年发现该方法是基于变分原理的里兹法的另一种形式，确立其理论基础。我国冯康在同一时期独立提出并证明了该方法：Melosh证明了位移法就是基于最小势能原理的Rayleigh-Ritz法冯康，基于变分原理的差分格式，应用数学和计算数学，1965，2(4):238-2621960至今：实际工程应用：平面空间板壳；静力动力、波动稳定；弹性塑性粘弹性、复合材料；固体流体、传热等连续介质力学；计算分析优化设计、与CAD技术结合。E.L.Wilson：编写了第一个公开的有限元软件SAP；通用有限元软件：SAP、ADINA、NASTRAN、ANSYS、ABAQUS、MIDAS等从半个多世纪以来有限单元法的萌芽、理论依据的证明和充实及其逐步的广泛应用可以看到，它的发展和计算机软硬件的发展基本上是同步的。如果没有计算机的强大软硬件支撑，有限单元法只有其微不足道的一点理论上的意义，而没有更为重要的实际应用的意义。2.1.2有限单元法概念(1)离散化离散化的过程是将连续体划分为有限数目、有限大小的单元的集合体。单元与单元之间只在指定点（即结点）连接，其他位置则一般保持连续即可。单元可以具有不同的形状，即单元外形可以不同；单元与单元之间可以有不同的连接方式，即单元的结点数目、位置可以不同。连续体连续体g图2.1连续体离散为单元集合体示例(2)单元分析对典型单元假设位移模式（由各结点位移插值），再分析单元的力学特性，建立单元的结点力与结点位移之间的关系，即单元刚度方程：{F}e=[k]{}e）并将各类荷载变换为作用在结点上的等效结点荷载。(3)整体分析将各单元刚度方程集成整体结构的整体刚度方程：{F}=[K]{}根据结点的平衡条件，得最终的有限元方程：[K]{}={R}求解该方程可得到未知的结点位移。(4)再次单元分析求出各单元的应变和应力。2.2弹性力学平面问题的矩阵描述2.2.1两类平面问题弹性力学的平面问题可分为平面应力问题和平面应变问题两类。实际上，所有的弹性力学问题都是空间问题。所谓平面问题，并不是说这个问题所分析的对象本身（如形状、荷载分布）是平面的，而是指该问题的形状、外部作用以及问题的解答（即由此产生的效应，如位移、应力等）只在平面内有变化，而沿着平面外就保持不变了。因此可以肯定地说，所谓的平面问题就是一个特殊的空间问题。那么，是不是一个问题的形状和外部作用（即已知的位移和应力边界条件）只在平面内发生变化，而沿着平面外保持不变了，这个问题就是平面问题呢？不是的，还必须附加其他条件，这一结论才能成立。这个附加条件就是该问题沿平面外的尺寸与平面内尺寸相比要么非常小（如无限短），要么非常大（如无限长）。如果符合前者条件，则弹性体内只存在平面内的应力，而平面外的应力均为零，故这类问题称为平面应力问题；如果符合后者条件，则弹性体内只存在平面内的位移或平面内的应变，而平面外的位移及应变均等于零，故这类问题称为平面应变问题。zzxt/2Oyt/2y体积力表面力图2.2平面应力问题示例xxOyz图2.3平面应变问题示例2.2.2(1)基本量应力：{}=[xyxy]T，应变：{}=[xyxy]T，位移：{f}=[uv]T；体积力：{p}=[pxpy]T，表面力：{q}=[qxqy]T。(2)基本方程几何方程：物理方程：{}=[D]{}，式中[D]：弹性矩阵。对平面应力问题，平衡方程的弱形式——能量原理1)虚功原理：弹性体在外力作用下处于平衡状态的充分和必要条件是，外力在满足变形协调推进的虚位移上所做的虚功等于其内力在相应的虚应变上所做的虚功，即2)最小势能原理：满足几何条件的各种可能位移中，真实位移使取驻值，即=(U+UR(3)平面问题的常用单元：三结点三角形单元六结点三角形单元四结点矩形单元四结点任意四边形单元八结点曲边四边形单元图2.4平面问题的常用单元单元的加密方法可分为p型加密和h型加密。前者是保持单元的大小和形状不变，而提高单元的插值函数的阶数；后者是指采用同一类单元，但加密单元的网格，即减小单元的尺寸，增加离散单元的数目。2.3三结点三角形单元分析2.3.1离散将一平面结构离散为ne个单元，设结点综述为n个，则结构的基本未知量取为这n个结点的位移，即{}=[u1v1u2v2…unvn]Tyyui(Ui)ixvi(Vi)uj(Uj)jvj(Vj)um(Um)mvm(Vm)图2.5三结点三角形单元2.3.2单元位移模式单元结点力向量：{F}e=[UiViUjVjUmVm]T单元结点位移向量：{}e=[uiviujvjumvm]T(1)位移模式的建立所谓位移函数或称位移模式，是指利用单元的结点位移将单元内任一点的位移用坐标的函数表示出来。对于三结点三角形单元，假设其位移模式是坐标的线性函数，则有u=1+2x+3yv=4+5x+6y系数1~6由6个结点位移分量（ui、vi、uj、vj、um、vm）确定。将位移模式写成结点位移的显式：u=Niui+Njuj+Nmumv=Nivi+Njvj+NmvmNi、Nj、Nm称为形状函数（形函数）或插值函数，其中ai、bi、ci是分母行列式第1行各元素的代数余子式，展开后可表示为ai=xjymxmyj，bi=yjym，ci=xj+xm（i,j,m轮换）jjyimxPj1imNi(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}=[N]{}e代入几何方程，得单元应变为：{}=[B]{}e[B]称为单元应变矩阵，[B]=[BiBjBm]，而各子块为。可见[B]为一常量矩阵，故三结点三角形单元为一常应变单元。将{}=[B]{}e代入物理方程，可得{}=[S]{}e=[D][B]{}e[S]称为单元应力矩阵，[S]=[S1S2S3]，对平面应力问题。可见[S]也是一常量矩阵，故三结点三角形单元为一常应力单元。2.3.4有限元方程的建立将连续体近似看作为由一系列只在结点相连的单元组装而成的集合体，并且对各个单元均建立了其位移模式，那么该集合体的最终位移法方程，即最终的有限元方程可由变分原理（此时为最小势能原理）建立。位移模式确定后，离散结构的可能位移就是由不同结点位移决定的分片函数。各可能位移函数中，真实的位移函数使离散结构的总势能取驻值（处于稳定平衡状态时为最小值），即=(U+UR令：{}e=[G]{}这里[G]是一62n的位置矩阵，表示单元结点位移{}e在整体结点位移向量{}中的位置；[k]e为单元刚度矩阵；[R]e单元等效结点荷载向量。于是有令：这里[K]为整体刚度矩阵；[R]为整体等效结点荷载向量，则有：由势能驻值原理=0可得或2.3.5单元刚度矩阵(1)单刚列式其中，(2)单刚性质单元刚度元素kpq的物理意义为：第q个结点位移分量为单位位移（其它结点位移等于0），所引起的第p个结点力分量。例如要获得元素k26的直观解释，可先将单元的所有6个结点位移约束住，即在每个结点处沿两个方向同时施加链杆约束，如图2.7所示。ymxk26iymxk26ij1单元刚度矩阵[k]的性质：1)对称性：kpq=kqp；2)奇异性：因具有刚体位移；3)每行（列）元素之和为零；例如，第6列元素的意义：当第六个结点位移等于1（6=1），其他结点位移等于0时，在各结点位移方向施加的结点力的大小，即图2.7中6个附加链杆约束上的约束力。因单元在这些结点力作用下处于平衡，故有：k16+k36+k56=0；k26+k46+k66=0。4)主对角元素恒大于零：kqp>0；5)[k]取决于单元的形状、方位和弹性常数，与所在位置（即平移或n转动）无关。(3)推导单刚的另一种方法——由单元的平衡条件导出物理方程物理方程几何方程位移模式虚功方程{f}=[N]{}e{}=[B]{}e{}=[S]{}e，[S]=[D][B]{F}e=[k]{}e，[k]=[B]T[D][B]tA结点位移位移应变应力结点力{}e{f}{}{}{F}e图2.8单元分析中各物理量的联系图假设单元发生了虚位移，则结点力在结点虚位移上做的虚功=单元应力在虚应变上的虚功（虚变形能）：令：又因{}e是任意的，故有{F}e=[k]e{}e2.3.6单元等效结点荷载在将连续体用近似的单元集合体代替以后，单元集合体中的每个结点可认为一方面受到外荷载的作用，另一方面受到连接该结点的各个单元的对其产生的约束力（即结点力）的作用。各结点在这两组力的作用下总是处于平衡的。但是，如果外荷载不是直接作用在结点上，而是作用在单元内部，那么怎么建立结点的平衡关系呢？一种有效方法是先建立单元的等效结点荷载，再对结点建立平衡关系。yYiyYiixXiYjjXjYmmXmPxPy图2.9单元的等效结点荷载分量2.3.7整体刚度矩阵(1)建立有限元方程的方法：1)由最小势能原理建立：[K]=([G]e)T[k]e[G]e2)由结点平衡条件建立(2)整体刚度矩阵[K]的集成方法：对号入座。[K]=([G]e)T[k]e[G]e，其中([G]e)T决定[kpq]e在[K]中的行，[G]e决定在[K]中的列位置。(3)整体刚度矩阵的性质[K]元素的物理意义Kpq：结构第q个结点位移分量为单位位移（q=1，其它结点位移=0），在第p个结点位移方向所施加的结点力的大小。如K25为结构第5个结点位移分量为单位位移，即5=1，其它结点位移均为零时，在第2个结点位移方向所需施加的结点力的大小。[K]的性质：1)对称性：Kpq=Kqp；主对角元恒大于零；2)稀疏矩阵，且一般呈带状分布；例如，平面问题最大半带宽=2(单元结点号之差最大值+1)。3)引入位移边界条件前为奇异矩阵，引入后为正定矩阵。2.3.8位移边界条件的引入以上集成的总体整体刚度矩阵并不包含位移边界条件的信息，这种刚度矩阵常称为原始刚度矩阵，它是一个奇异矩阵。要使得最终的有限元方程可解，必须引入唯一边界条件，排除结构的刚体位移。引入位移边界条件的常用方法有以下几种：(1)直接引入法（矩阵缩小法）将[K]{}={R}改写为如下子块形式：其中{a}、{b}分别为未知结点位移向量和已知结点位移向量，[Kaa]、[Kab]、[Kba]、[Kbb]、{Ra}、{Rb}分别为与之对应的刚度矩阵和荷载向量子块。将上述方程的上半部分（对应已知位移）展开，得[Kaa]{a}={Ra}[Kab]{b}由该方程可解出未知结点位移。该方法将改变原始刚度矩阵的阶数及元素的顺序，不利于程序实现，因此在计算机编程中一般很少采用。(2)对角元素改1法（零位移边界）设结构的总自由度数（即结点位移总数）为N，若第i个结点位移为零（即i=0），则将[K]中对角元素Kii改为1，而第i行和第i列的其他元素改为0，荷载向量[R]中的第i个元素也改为0。其实质是将原有限方程中的第i个方程用方程i=0代替，而其他方程中与i对应的系数也改为0，表明i对其他方程没有影响，同时保证了修改后的刚度矩阵仍具有对称性。(3)乘大数法若第i个结点位移已知，即，则将整体刚度矩阵[K]中的对角元素Kii改为Kii，其中为一大数（如1020），而荷载向量[R]中的第i个元素Ri改为，原方程成为以下形式其实质是将原有限方程中的第i个方程用的以下近似方程代替：将上述方程各项同除以大数，除第i项及方程右端项外，其他各项均趋于0，故等价于。2.3.9位移模式与解答的收敛准则(1)有限元解答的收敛准则为使有限元解答能够收敛于精确解，单元位移模式应满足以下条件：1)位移模式必须包含单元的刚体位移；对弹性力学平面问题，其刚体位移表达式为：u=u0y，v=v0+x因此，平面问题的位移模式必要包含上述刚体位移表达式中的各项，才能保证最终解答的收敛性。2)位移模式必须包含单元的常量应变；条件(1)、(2)合并起来可称为完备性要求。对平面问题来说，就是要求位移模式必须包含常数项和完整的一次项。完备性条件是解答收敛的必要条件。3)位移模式应尽可能保证位移的连续性。该条件实际上就是协调性条件，但一般情况下并不是一个必要性条件。如果位移模式同时满足上述完备性和协调性条件，那么就组成了解答收敛的充分条件。对一般的弹性力学平面或空间问题，只需要求单元内部以及相邻单元的公共边界上的位移本身连续，即位移模式具有C0连续性。对有些问题，可以放松对协调性的要求，只要通过分片试验，那么也能保证解答的收敛性。三结点三角形单元的位移模式既满足完备性，又满足协调性的要求（在单元边界上呈线性分布，可由两个结点位移唯一确定），是一种协调单元。证明如下：单元内：单值连续；yx(1)ijmp(2)相邻单元之间：uij(1)=uij(2)？vijyx(1)ijmp(2)ij边的方程：y=ax+b，则uij=1+2x+3(ax+b)=cx+duij(1)、uij(2)均为坐标的线性函数，故可由i、j两点的结点位移唯一确定。图2.10两相邻单元(2)多项式位移模式的一般选择规则位移模式应与单元局部坐标的选取无关，即满足所谓的几何各向同性。对于平面问题，位移模式中的x和y的各阶项应保持对称，即有了xmyn项，则应同时具有xnym项。为保证位移模式关于x和y坐标的对称性，通常从以下的Pascal三角形中选取多项式位移模式的各项：1xyx2xyy2x3x2yxy2y3x4x3yx2y2xy3y4…………图2.11Pascal三角形根据最小势能原理建立的有限元，是以结点位移作为基本未知量，这种有限单元称为位移元。由位移元得到的近似解答总体上是精确解的一个下限。2.4高精度的三角形单元、矩形单元2.4.1高精度三角形单元(1)六结点三角形单元（二次单元T6）位移模式取坐标的完整二次式：u=1+2x+3y+4x2+5xy+6y2v=1+2x+3y+4x2+5xy+6y2该位移模式包含了常数项和完整的一次项，满足完备性要求；在单元的边界上位移呈二次抛物线分布，可由三个结点位移唯一确定，故又满足协调性的要求，是一种协调单元。(2)十结点三角形单元（三次单元T10）位移模式取坐标的完整三次式：u=1+2x+3y+4x2+5xy+6y2+7x3+8x2y+9xy2+10y3v=1+2x+3y+4x2+5xy+6y2+7x3+8x2y+9xy2+10y3该位移模式包含了坐标的完整一次式（常数项和纯一次项），满足完备性要求；在单元的边界上位移呈三次曲线分布，可由4个结点位移唯一确定，故又满足协调性的要求，是一种协调单元。ijijmijm3122.12高精度三角形单元(3)面积坐标表示的三角形单元形函数在推导三角形单元的列式时，直接利用整体坐标系（为直角坐标）下的位移模式将使得运算十分繁琐和复杂。如果采用三角形单元内的一种局部坐标——面积坐标作为自然坐标，则可以使列式推导大为简化。定义单元内任一点P的无量纲面积坐标（Li,Lj,Lm）为i(1,0,0)Pj(0,1,0)m(0,0,1)AiAjAmLii(1,0,0)Pj(0,1,0)m(0,0,1)AiAjAm各种单元的形函数：3结点三角形单元（线性单元T3）：Ni=Li（i,j,m）6结点三角形单元（二次单元T6）：Ni=(2Li1)Li（i,j,m）N1=4LjLm,（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=1+2x+3y+4xy1yx234ov=1+2x+1yx234o可保证位移模式的完备性和协调性。若写成形函数形式，则为2.14四结点矩形单元这里的Ni（i=1,2,3,4）可以先求出8个待定系数在获得，也可以根据形函数的性质直接构造，例如对N1，可设该表达式满足在结点2、3、4取值均为零的性质；再令Ni(1)=1，可得待定系数这里设矩形单元的边长各为2a、2b(2)局部坐标下的形函数设矩形长和宽各为2a和2b，在矩形形心为原点建立局部坐标（自然坐标）系o4123=1=1=1=12.15局部坐标系下的矩形单元则四个角点的局部坐标分别为1(1,1)，2(1,1)，3(1,1)，4(1,1)。位移模式为：Ni（i=1,2,3,4）可根据形函数的性质直接构造出来：或(3)单元应变矩阵和应力矩阵单元应变为：{}=[B]{}e其中[B]=[B1B2B3B4]，。单元应力为：{}=[S]{}e=[D][B]{}e其中[S]=[S1S2S3S4]，[Si]=[D][Bi]2.5C在有限单元法中，根据结点布置方式的不同可将矩形单元分为两类，一类称为Lagrange矩形单元，另一类称为Serendipity矩形单元。前者在单元纵横网格线的交点上均布置结点，而后者仅在单元的边界线上布置结点，如图2.17所示。(a)Lagrange矩形单元(b)Serendipity矩形单元2.16两类矩形单元2.5.1Lagrange矩形单元(1)一维Lagrange插值多项式过n个结点（坐标分别为x1,x2,…,xn）的一维Lagrange插值多项式为xxx1x2x3xnxi12.17一维Lagrange插值多项式l1(x)例如n=2，则有：。若令x1=0，x2=l，则。用该多项式作为形函数，可满足形函数的性质。(2)Lagrange矩形单元的形函数在矩形的各个网格交点上均布置结点，如水平方向r+1个，竖向p+1个。于是第I列J行结点i的相应形函数：Ni=NIJ=lI(r)()lJ(p)()其中，lI(r)()、lJ(p)()均为一维Lagrange插值多项式，Ni在结点i为1，在其他结点处为零；单元边界上的结点数=形函数阶数+1，能够保证边界位移的唯一性和协调性。这类单元包含较多的内部结点，增加了单元的自由度，而实践证明这些自由度的增加通常并不能有效提高单元的精度。对nn结点的Lagrange单元，其Pascal三角形中包含的项数如图2.19所示。(I(I,J)(r,p)(0,0)11lI(x)lJ()1xnynxnynxx2x3yy2y3图2.18Lagrange矩形单元及插值模式图2.19Pascal三角形包含项数2.5.2.Serendipity矩形单元(1)结点布置特点这类单元只在其边界上布置结点，但不同边界上可布置有不同数目的结点。(2)形函数的构造方法如果只有四个角点上布置结点（R4单元），则如果增加结点5使之成为5结点矩形单元（或称R5单元），则满足(N5)j=5j（j=1,2,3,4,5）。原不满足，需对其进行修正才能满足：(N1)5=0，(N2)5=0，且保证N1，N2在除5结点外的其他4个结点取值不变。412351011011/21N5图2.20五结点矩形单元图2.21形函数N1的修正对于8结点Serendipity矩形单元（R8单元），利用上述方法，可得到前4个结点对应的形函数分别为或写成统一表达式：后4个形函数分别为对于p次单元的边界内结点，其相应的形函数为（或）的p次和（或）的一次Lagrange多项式的乘积；而角结点的形函数为四结点单元的相应函数与各内结点形函数乘以一分数的差。利用该方法可以很方便地构造一些过渡单元的形函数。2.6平面等参单元直接用形状规则的单元对复杂形状的结构进行离散绘遇到边界难于处理的问题，因此应设法采用适当的方法对规则单元进行坐标变换，使之转化为斜边或曲边的单元。有限元中最常用的变换方法是等参变换，即坐标变换采用与单元位移插值一致的形函数。通过等参变换的单元称为等参单元。2.6.1四结点四边形等参元对于4结点任意四边形单元，若采用双线性的位移模式，则位移在单元斜边界上呈二次抛物线分布，由两个结点位移不能唯一确定。为此在单元内设法建立一个局部坐标（自然坐标），使得单元各边界的局部坐标分别为一常量（或1），这样，如果位移模式采用局部坐标的双线性函数，则能够满足公共边界位移的唯一性。如果单元边界的局部坐标分别为1，则在局部坐标下，单元的形态就是图示的4结点正方形单元，该单元称为基本单元或母单元，其形函数为。位移模式为：1y1yx2344123=1=1=1=1(a)基本单元(b)实际单元2.22四结点等参数单元如果同时利用上述形函数作为局部（自然）坐标(,)向整体（直角）坐标(x,y)的变换函数，则可以建立两种坐标之间的映射关系如下：根据形函数的性质：(Ni)j=ij，Ni=1，容易验证上述变换式在四个结点处给出了节点的整体坐标，而在单元的四边上，其中一个局部坐标为1，另一局部坐标按线性变化，因而正好给出了整体坐标下四边的直线方程（由两个结点整体坐标唯一确定）。例如对12边，因此上述变换式正确地反映了局部（自然）与整体（直角）坐标之间的映射关系。经变换后的实际单元称为子单元。子单元的位移模式仍为：利用等参变换可以构造8结点、12结点、20结点等更高次的四边形曲边等参单元（参见图2.23）。o15247386yx15247386图2.238结点曲边四边形等参单元2.6.2局部与整体坐标的微分和积分变换式根据复合函数的求导法则，写成矩阵形式，式中J称为雅可比（Jacobi）矩阵，对于具有m个结点的平面等参单元，若反过来用局部坐标表示整体坐标，则可对上式作反变换，式中，Jij*是J各元素Jji的代数余子式；J是J的行列式，称为雅可比（Jacobi）行列式。面积微分的变换：dA=dxdy=Jdd。2.6.3单元刚度矩阵和单元等效结点荷载向量对于一具有m个结点的平面等参单元，其应变向量可写为{}=[B]{}e其中，单元应变矩阵：[B]=[B1B2…Bm]。单元应力向量：{}=[S]{}e=[D][B]{}e单元刚度矩阵：单元等效结点荷载列阵：式中{Rp}e、{Rq}e分别为体积力和表面力引起的等效结点荷载，且对其中一个局部坐标（或）为常数的边界，线积分ds为或上述积分式通常采用Newton-Cotes或高斯（Gauss）数值积分求得。可见，等参单元的刚度矩阵、等效结点荷载列阵等的计算都在规则的母单元中进行，因此各种形式的积分运算都可以采用标准的数值积分方法进行，使得不同工程问题的有限元分析能够采用统一的通用化程序实现。2.6.4等参变换的应用条件等参单元的应用条件是母单元与实际单元之间存在一一对应关系。具体到局部与整体坐标的变换式，要求变换矩阵即雅可比[J]非奇异，雅可比（Jacobi）行列式J0。为确保坐标变换一一对应，在单元划分时应避免：使任意两个结点退化为一个结点而使d=0或d=0；还应避免因单元过于歪斜而导致d与d发生共线。123=1=1=1=14123,4=1=1=1=1(a)结点退化(b)单元边界接近共线2.24畸形等参数单元2.6.5例题分析例2.1图示悬臂平面结构，长2m，高1m，试用图示单元进行分析。112346(1)(2)(3)5(4)xyPjjmi(1)(3)jmi(2)(4)(a)基本单元(b)实际单元2.25例2.1图（三角形单元）1.三结点三角形单元分析(1)离散化：划分为4个单元，共6个结点。(2)依据图示的单元局部编号规则，4个单元的刚度矩阵均相同，为为简便起见，设泊松比=0，于是有(3)根据局部与整体结点编号的对应关系集成整体刚度矩阵：(4)引进位移边界条件2.四结点四边角形等参单元分析112346(1)(2)5xyP41232.26例2.1图（四边形单元）计算步骤：(1)离散化：划分为2个单元，共6个结点。单元上下短边0.75m，长边1.25m。(2)求，[J]，|J|，[J]-1；(3)求[B]在各积分点的数值[Bi]g；(4)利用高斯积分计算并形成[k]；(5)集成[K]、[P]；(6)引进位移边界条件。112P例2.2牛腿受竖向集中力作用，且结点1、2发生已知的水平位移，求应力。2.27例2.2图2.7有限元程序实现有限元程序——前处理程序：生成网格及数据文件主体分析程序：核心计算分析后处理程序：进行结果处理，生成可直接使用或查看的结果文件、图形文件。2.7.1有限元程序设计的一般步骤算法描述和列式推导；框图设计；代码编写；上机调试、考核；编写应用说明；修改、补充、完善。程序设计的一般要求：具备较齐全的功能、较强的通用性和可移植性、较好的可扩充性、良好的可读性、足够的可靠性、良好的自适应性。2.7.2输入数据及分类1.控制数据：结点总数、单元总数、约束总数、荷载总数、问题类型数等2.几何数据：结点坐标、单元信息（各单元的结点编号）、约束条件、单元类型数（弹性模量E、泊松比、厚度t不同为一类）3.物性数据：弹性模量E、泊松比、厚度t4.荷载数据：荷载类型（集中、分布）、位置、方向、大小等2.7.3三结点三角形单元分析平面问题的主体程序1.程序框图（参见图2.28）2.结构化程序设计方法模块化——由1个主程序和若干子程序组成。子程序由通用性，采用可调数组；主程序采用动态数组存储技术3.动态数组存储技术（1）按整型和实型定义两个大型共享数组，如A(1000000)、M(1000000)；（2）设计动态数组表：将各子程序中的变界（可调）数组按各自实际需要的大小分配一维数组的空间；（3）动态数组覆盖技术：全部或部分覆盖。读入控制数据读入控制数据开始读入几何、物性、荷载数据平面应力/应变?E=E0/(1-02),=0/(1-0)E=E0,=0计算单元刚度元素叠加到整体刚度矩阵中计算等效结点荷载引入位移约束条件解线性代数方程得结点位移计算单元应力、反力等输出结果结束图2.28平面问题主体分析程序框图2.8平面杆系结构的有限元2.8(1)基本方程xxyp(x)q(x)图2.29受任意荷载的等截面直梁几何关系：写成矩阵系数，，{f}=[uv]T，内力-位移关系：平衡关系：边界条件：(2)离散将一平面杆件结构离散为ne个单元，n个结点，基本未知量为结点位移{}=[u1v11u2v22…unvnn]T对应的结点力向量为{F}=[X1Y1M1X2Y2M2…XnYn(3)位移模式局部坐标下的单元结点位移向量和单元结点力向量：xxyx’y’viuivjujij图2.30等截面直梁单元由梁的平衡方程可知，在结点荷载作用下梁的轴向位移沿梁轴呈线性分布（因只有杆端作用轴向荷载），横向位移呈三次曲线分布（因弯矩沿杆轴成线性分布），故假设u=1+2xv=3+4x+5x2+6x3将结点i,j的位移代入可求出6个待定系数1~6。将单元位移写成结点位移的显式，有u=N1ui+N4ujv=N2vi+N3i+N5vj+N6j位移模式的矩阵表示u、v独立插值，但不独立插值，故要求C1连续：不仅u、v本身连续，v的一阶导数也要连续。(4)单元应变矩阵和应力矩阵将位移模式{f}=[N]{}e代入几何方程，得单元应变为：将{}=[B]{}e代入物理方程，可得{}=[D][B]{}e,其中(5)单元刚度矩阵局部坐标下的单元刚度矩阵：单元刚度方程为：{F}e=[k]{}e[k]中的每一元素krs称为单元刚度系数，其物理意义是：第s个杆端位移为单位位移，而其他杆端位移为零时所引起的第r个杆端力。单元刚度矩阵是一个对称矩阵，即krs=ksr，这可由反力互等定理证明；单元刚度矩阵还是一个奇异矩阵，这是由于单元中包含刚体位移。局部向整体坐标的变换：{’}e=[T]T{}e{F’}e=[T]T{F}e[k’]=[T]T[k][T]式中，[T]为一66的坐标变换矩阵。假设自整体坐标Ox’轴沿转角正方向（即顺时针方向）转到局部坐标Ox轴的角度为，则{F’}=[k’]{’}e坐标变换矩阵[T]是一个正交矩阵，即[T]1=[T]T，故结点位移和结点力由整体向局部坐标的变换式为{}e=[T]{’}e，{F}e=[T]{F’}e

2.8基本假定：垂直于梁轴线的截面在变形后仍保持为平面，但不再垂直于变形后的轴线。挠度由两部分组成：弯曲变形部分和剪切变形部分，即v=vb+vs；dxdx杆轴线斜率：dv/dx=+（dxdx图2.31发生弯曲和剪切变形的直梁微段(1)在经典梁单元基础上引入剪切变形的影响vb用三次式插值，参见上一节。vs用线性插值：vs=N7vsi+N8vsj。由此可导得单元刚度矩阵为式中，为剪切影响系数，k为截面剪应力不均匀分布修正系数，对矩形截面k=1.2，圆形截面k=10/9。(2)Timoshenko梁单元对u、v、均独立插值：u=N1ui+N2u2，v=N1v1+N2v2，=N11+N22式中：截面的转角，dv/dx：杆轴线的斜率。当l/h时只能得到零解，即将出现剪切自锁（Shearlocking）。这是由于v、同阶插值（实际上放大了剪切应变项的量级），故dv/dx与不同阶，从而导致剪应变=dv/dx不能处处满足（dv/dx为常数，为一次式），除非=常数，意味着梁不能发生弯曲。解决方法：减缩积分——数值积分时采用比精确积分要求少的积分点数，例如对两结点梁单元采用一点积分；减缩积分后[ks]中的l2/3和l2/6项均改为l2/4。假设剪切应变——对剪应变另行假定插值形式；替代插值函数——计算剪应变时，对采用低一阶的插值函数，如两结点梁单元其插值函数为常数1/2，即。两结点Timoshenko梁单元包含横向刚体位移（v=c）和刚体转动（=dv/dx=c），包含常剪切应变（=0，dv/dx==c），但不包含常弯曲应变的位移状态（=cx，v=0.5cx2），不能分析纯弯问题（将伴随剪切应变），因而是一种较低级的C0连续型梁单元，实际计算中多采用三结点或更多结点的梁单元。2.9板弯曲的有限元2.9.1薄板小(1)基本假定1）直法线假定：z=0，zx=0，yz=0；2）不计z引起的变形（物理方程与平面应力问题相同）；3）小挠度假定：变形后中面无面内位移。基本方程几何方程：{}=[L]w，，式中{}为曲（扭）率向量。应变：{}=[xyxy]T=z{}。物理方程：{M}=[MxMyMxy]T=[D]{}式中[D]为板弯曲的弹性关系矩阵，,为板的弯曲刚度。平衡方程：或边界条件：2.9(1)单元的自由度单元与单元之间只在结点连接，由于相邻单元之间要传递横向力及力矩，故结点须作为刚接。每个结点有3个自由度：挠度w、角位移x、y。结点i的位移：（i=1,2,3,4）相应结点力：（i=1,2,3,4）单元结点位移向量：{}e=[1234]T边界条件如下：结点i双向铰接——wi=0结点i双向刚接——wi=xi=yi=0绕x（或y）转向铰接，绕y（或x）转向刚接——wi=yi=0（或wi=xi=0）xzxzy=xyzw1x1y11234图2.32矩形博班单元(2)单元位移模式w1因薄板的位移、变形及内力等均可用挠度w表示，故只需对w建立位移模式。w1w=1+2x+3y+4x2+5xy+6y2+7x3+8x2y+9xy2+10y3+11x3y+12xy3将4个结点的位移代入可求出12个待定系数1~12。将单元位移写成结点位移的显式，有w=[N]{}e，[N]=[N1N2N3N4]，设单元在x、y方向的边长分别为2a和2b，则在局部坐标系o其中局部与整体坐标之间的关系为这里xc、yc为单元中心的x、y坐标。(3)位移模式的完备性和协调性分析完备性分析：前三项1+2x+3y代表了薄板单元的刚体位移（1代表薄板沿z方向的刚体平移，2、3分别代表薄板绕y和x轴的刚体转动）；4x2+5xy+6y2代表单元的常应变（常曲率和常扭率）。连续性分析：要求C1连续。例如在y=常量的单元边界12上，w为x的三次式，即：w12=1+2x+3x2+4x3，可由w1、w2、y1、y2唯一确定。该边上x也是x的三次式（该三次式与w的三次式相互独立），即(x)12=(w/y)12=5+6x+7x2+8x3，但只有x1和x2为限制条件，故不能唯一确定，为一非协调元。该单元经验证能够通过分片试验，即当单元趋于无限小时能保证应变为常量，解答收敛于正确解。分片试验：任取一至少包含一个内结点的单元片，赋予各结点与常应变状态相应的位移值，验证由此而产生的结点力是否处于平衡，若是则认为单元片通过分片试验。2.9位移模式：w=1+2x+3y+4x2+5xy+6y2+7x3+8x2y+9xy2+10y3因只有9个自由度，故上式需删去一项，如将8、9两项合并或舍去xy项，均不能保证收敛，因此改用面积坐标表示。面积坐标的性质：1=Li+Lj+Lmx=xiLi+xjLj+xmLmy=yiLi+yjLj+ymLm(a)x、y一次式：Li、Lj、Lm(b)x、y二次式：Li2、Lj2、Lm2、LiLj、LjLm、LmLi(c)x、y三次式：Li3、Lj3、Lm3、Li2Lj、Lj2Lm、Lm2Li、LiLj2、LjLm2、LmLi2、LiLjLmx、y的完全三次式应至少包含（c）中的4项及（a）、（b）、（c）中的6项共10项的线性组合。连续性：单元边界上，w为三次式，可由两端点的w及w/s唯一确定，故能够保证连续；w/n为二次式，由两端点的w/n不能唯一确定，故w/n不能保证连续。2.9构造方法：1.保持结点参数不变，采用附加校正函数法使单元边界协调2.增加结点参数，例如包含w的二阶或更高阶导数，如21自由度的三角形单元3.杂交混合元、拟协调元、广义协调元等。2.9(1)Mindlin-Reissner中厚板理论1）z=0，zx0，yz0（放弃了直法线假定）；2）不计z引起的变形（物理方程与平面应力问题相同）；3）变形后中面无面内位移（小挠度假定）。(2)四结点矩形单元对挠度w和转角x、y均独立插值，其位移模式为：该单元只要求C0连续。当用于薄板问题时会出现剪切自锁（Shearlocking），这是由于不满足Kirchhoff假定而出现的剪应变能失真。解决方法：减缩积分；假设剪切应变；替代插值函数。

习题图示为一平面应力状态的三结点直角三角形单元，厚度t，弹性模量E，剪切模量G=E/[2(1+)]，设泊松比=0，结点坐标如图。若采用线性位移模式（位移函数），试求出：xym(0,1)jxym(0,1)j(1,0)i(0,0)(2)应变矩阵[B]；(3)应力矩阵[S]；(4)单元刚度矩阵[k]；(5)[k]的每行之和及每列之和，并说明其物理意义。题2.1图为使有限单元解收敛于正确解，位移模式应满足那些条件？对于平面四结点矩形单元，若位移模式取为：u=a1+a2x+a3y+a4x2，v=b1+b2x+b3y+b4y2，试分析该位移模式是否满足这些条件，并说明具体理由。为使有限单元解收敛于正确解，位移模式应满足那些条件？四结点矩形薄板单元具有12个自由度，其位移模式取为：w(x,y)=1+2x+3y+4x2+5xy+6y2+7x3+8x2y+9xy2+10y3+11x3y+12xy3，试分析该位移模式是否满足这些条件，并说明具体理由。形函数有哪些主要性质？试由这些性质直接构造图示六结点矩形单元的形函数，写出单元中心点P(a/2,b)处的位移用结点位移表示的表达式。114bxyba2365mmxyji题2.4图题2.5图图示为平面问题的一个三结点三角形单元。(1)试问单元刚度矩阵[k]有哪些主要特性？其依据各是什么？(2)附图说明[k]元素k52的物理意义。(3)[k]的每行之和及每列之和各为何值，其物理意义是什么？图(a)所示的平面连续体结构已划分为两个三角形单元，在图(a)坐标系及图(b)局部编号下，两单元的刚度矩阵左下子块均为：。附图说明单元(1)的刚度元素k36的物理意义；试由上述单元刚度矩阵子块形成结构的总体刚度矩阵；分别采用手算方法和一种计算机方法引进图中的位移边界条件，写出图示荷载作用下的最终有限元方程；假设结点位移v2、u3、v3、u4均已求得(作为已知)，试在此基础上求出结点2和结点4的支座反力。44123(1)(2)Pxyjjmiij(1)(2)m(a)(b)题2.6图利用形函数的性质，直接构造出图示六结点正方形单元（边长为2）的形函数，写出单元中心点o的位移用结点位移表示的表达式mxmxyjiPyPx1111152463题2.7图题2.9图有限单元法中，一个二维单元在坐标平面内分别发生平移和转动，单元刚度矩阵[k]是否发生改变？为什么？应力矩阵[S]又如何变化？在有限元分析中，非结点荷载需移置为等效结点荷载，移置的原则是什么？试根据该原则，导出三结点三角形单元内任一点（x，y）处作用集中荷载{P}=[Px,Py]T时的等效结点荷载表达式。已知形函数矩阵为[N]，结点位移向量为{e}。三结点三角形单元的材料容重为，厚度为t，试导出单元在自重作用下的等效结点荷载向量。三结点三角形单元的ij边作用一法向的线性分布荷载，其中i、j点的集度各为qi,qj，试导出单元的等效结点荷载表达式。图(a)示平面结构已划分为两个四边形单元，各结点的整体编号如图所示。已知两单元的上下短边0.75m，长边1.25m，竖向垂直尺寸1m，其所对应的等参基本单元如图(b)所示。试求（=0,=0）点对应的单元(1)、(2)的整体坐标（x,y）已知各结点的位移ui,vi（i=1,2,3,4），均为已知，试求结构在（x=0.625,y=0）处的位移1(1,1)2(1,1)3(1,1)4(1,1)O(0,0)12346(1)(2)5yx(a)实际单元(b)基本单元题2.12图有限单元法中，引入位移边界条件前后结构的总体刚度矩阵各有哪些基本性质？引入位移边界条件的方法有那些？总体刚度矩阵在计算机中常用那些方法进行存储？各有什么优缺点？图示平面结构已划分为三结点三角形单元，试完成以下分析：怎样对结点和单元进行编号可使总体刚度矩阵的存储量最小？此时存储阶数为几乘几？若在结点C处添加一刚度系数为k0的水平弹簧支座（弹簧不作为单元），试根据总体刚度元素的物理意义说明此时总刚矩阵有何变化？CCAB题2.14图在板弯曲问题的有限元分析中，有两类常用的板单元：基于经典薄板弯曲理论的板单元和基于中厚板（Mindlin板）理论的板单元。这两类单元的基本假定、连续性要求有何不同点？Mindlin板单元在用于薄板时会遇到什么困难？常用那些方法克服该困难？

第3章加权残数法3.1基本原理及发展概况3.1.1基本原理加权残数（余量）法直接从微分方程和边界（初值）条件出发，先假设一包含待定系数的近似试函数，将试函数代入微分方程和边界（初值）条件中，一般不能精确满足而出现误差值，称为残值（或余量）；使残值在计算区域和边界上按某种加权平均为零，从而确定出待定系数，得到最终的近似解。基本方程：L(u)p=0，V（域内）边界条件：B(u)g=0，S（边界）L、B为算子，p、g为已知函数。为求近似解，设试探函数代入方程和边界条件，得残值：RL=L(u)-p（域内），RB=B(u)-g（边界）令残值在某种平均意义（加权积分）上等于零，则有Wj、Wjs分别为域内残值和边界残值的权函数。由上式得到关于待定系数i的代数方程组，解得待定系数。3.2.2发展概况加权残数法的基本思想在19世纪初就已提出，到了1920年代用于求解微分方程，1950年代CrandallS.H.统一命名为“加权残数（余量）法”（WeightedResidualsMethod）。1970年代前主要应用于计算流体力学、空气动力学、热传导等领域，而较少用于固体力学领域。国内，徐次达教授于1970年代末1980年代初首先引进并应用于固体力学计算中。在1980-90年代得到了很大发展，自1982年起召开了多届专门的学术会议，在国内外期刊发表了大量的论文。涉及多方面的研究领域：静力、动力、稳定，线性、非线性，还于其他一些数值方法（如有限元法）、半解析法相互结合，用于求解复杂的工程问题，如加权参数法与有限元的耦合求解法、摄动加权参数法等。3.2加权残数法的分类3.2.1按试函数分类按试函数是否满足微分方程和边界条件分为：内部法、边界法和混合法。(1)内部法试函数满足边界条件而不满足微分方程，RB=0，(2)边界法试函数满足域内微分方程而不满足边界条件，RL=0，(3)混合法试函数既不满足域内微分方程而不满足边界条件，RL0，RL0，分析实际工程问题时，应尽可能采用内部法或边界法，使问题得到简化。最常用的是内部法（因边界条件相对易于满足），对较复杂的问题才用混合法。3.2.2按权函数分类按权函数分类，加权残数法可分为5种基本方法：子域法、配点法、最小二乘法、矩量法和Galerkin（伽辽金）法，现以内部法为例分别说明。(1)子域法（sub-domainmethod）在求解域V中取m个子区域Vj（j=1,2,…,m），使微分方程的残值在各子域内积分为零，即（试函数可包含m个待定系数）加权积分式：权函数：例3.1：微分方程，边界条件，求u的近似表达式。采用内部法。选取满足边界条件的试函数：u=x(1x1+2x+3x2+…)取前两项：u=x(1x1+2x)得残值：R=L(u)px+(2+xx2)1+(26xx2x3)2取两个子域，即m=2。子域1：x=0~0.25，有子域2：x=0.25~0.5，有解得：1=0.194132，2=0.191846与精确解的误差为3~4%。(2)配点法（collocationmethod）在求解域V中取m个点xj（j=1,2,…,m），使各选点的残值为零，即Rx=xj=0（j=1,2,…,m）得到关于m个待定系数的线性代数方程[C]{}={b}。该方法的权函数：Wj=(xxj)式中的(xxj)称为Diracdelta函数，简称函数。该函数的定义：，函数的性质：xx(x-)图3.1函数例3.2：用配点法求解例3.1所述问题。选取满足边界条件的试函数:u=x(1x1+2x)选点：x1=0.25，x2=0.5，得解得：1=0.193543，2=0.184332。与精确解的误差为2~3.6%。(3)最小二乘法（leastsquaresmethod）基本思想：选择试函数，使残值的平方和最小。有两种形式：连续型、离散型。(1)连续型建立域内残值平方的积分式：。为使I最小，取极值条件将前一式代入，有该式提供了关于m个待定系数的线性代数方程[C]{}={b}最小二乘法的权函数为：例3.3用最小二乘法求解例3.1所述问题。选取满足边界条件的试函数（两项近似）：u=x(1x1+2x)得残值：R=L(u)px+(2+xx2)1+(26xx2x3)2，解得：1=0.187542，2=0.169470与精确解的误差为1.5~2.5%。(2)离散型取n个配点，使各配点残值的平方和最小。其中的Ri（i=1,2,…,n）是j（j=1,2,…,m）的线性表达式。由极值条件可得：，经整理可得到关于{}的代数方程[C]{}={b}例3.4：用离散型最小二乘法计算例3.1。试函数选取及残值计算同前，即：u=x(1x1+2x)，得残值：R=L(u)px+(2+xx2)1+(26xx2x3)2取配点数n=3，并选x1=0.25，x2=0.5，x3=0.75，代入上述残值表达式可得R1、R2、R3以及等的表达式，代入离散型最小二乘法的方程式中，有解得：1=0.192973，2=0.172043，误差0.0~0.53%。如选配点x1=0.3，x2=0.6，x3=0.9，则误差0.72~0.97%。(4)矩量法（methodofmoment）权函数为Wj=xj（j=0,1,2,…），则有理论依据：有限矩量定理——若在有限区间[a,b]上连续的函数f(x)满足的n个关系式，则f(x)在[a,b]上至少变号n次，即f(x)=0在[a,b]上至少有n个相异的根，当n时，f(x)0。ff(x)xo图3.2变号n次的连续函数例3.5用矩量法计算例3.1。选取满足边界条件的试函数（两项近似）：u=x(1x1+2x)残值：R=L(u)px+(2+xx2)1+(26xx2x3)2取权函数为1，x，解得：1=0.187982，2=0.169492与精确解的误差为1.6~2.3%。(5)伽辽金法（Galerkinmethod）该方法1915年由俄国工程师Galerkin提出，是Ritz变分法的推广。取满足边界条件的试函数其中k（k=1,2,…,m）为基函数。取权函数为试函数的各项基函数，即Wj=j（j=1,2,…,m），则物理意义：使残值R与试函数的各项（即基函数）均正交。例3.6用伽辽金法计算例3.1。选取满足边界条件的试函数（两项近似）：u=x(1x1+2x)残值：R=L(u)px+(2+xx2)1+(26xx2x3)2取权函数W1=x(1x)，W2=x2(1x)，解得：1=0.192412，2=0.170732与精确解的误差为0.5%。5种基本方法可以混合使用，产生新的混合加权残数法，如最小二乘配点法（即离散型最小二乘法）、Galerkin配点法等。3.3伽辽金法的应用3.3.1Galerkin法解一维问题例如，图3.3所示简支梁受均布荷载，求梁的弹性挠曲线。微分方程：EIv(4)q=0，边界条件：(v)x=0或x=l=0，(v’’)x=0或x=l=0。yyxEI=常数lqAB图3.3简支梁受均布荷载取满足边界条件的试函数：，残值：权函数：Galerkin方程组：利用三角函数的正交性：完成上述积分计算，可得m个解耦的代数方程：解得：与精确解：相比，误差在0.5%以内。3.3.2矩形薄板弯曲的Galerkin解法周边固支的矩形薄板弯曲问题，当荷载情况较为复杂时，要获得问题的精确解比较困难，但利用Galerkin解法却十分方便。yyxBDAC22bO图3.4四边固支的矩形薄板微分方程：D4wq=0，边界条件：(w)x=a=0，(w)y=b=0，由对称性取满足边界条件的挠曲试函数为w(x,y)=(x2a2)2(y2b2)2(1+2x2+3y2+)取一阶近似：w(x,y)=1(x2a2)2(y2b2)2代入微分方程得残值R=8D1[3(y2b2)2+3(x2a2)2+4(3x2a2)(3y2b2)]q权函数W=(x2a2)2(y2b2)2代入Galerkin方程：对于均布荷载，可解得对于方板（a=b），中心挠度wmax=0.0213qa4/D，精确解wmax=0.0202qa4/D，误差5.4%。若取二阶近似，则误差可降低。3.4最小二乘法的应用3.4例.两端固定梁受均布荷载，求弹性挠曲线。微分方程：EIv(4)q=0，边界条件：(v)x=0或x=l=0，(v’)x=0或x=l=0。yyxlEI=常数qAB图3.5两端固定梁受均布荷载(1)内部法取满足边界条件的试函数：v=x2(lx)2残值：R=24EIq，权函数：W=R/=24EI加权残数方程：正好为一精确解，因试函数选择得恰当。(2)混合法取试函数为：v=0+1x+2x2+3x3+4x4既不满足方程，也不满足边界条件。取相对权函数W=1，则（j=0,1,2,3,4）内部残值：RL=24EIq，边界残值：RB1=v|x=00=0，RB2=v|x=l=0+1l+2l2+3l3+4RB3=v’|x=0=1，RB4=v’|x=l=1+22l+33l2+44l代入最小二乘法方程中，可解得故挠曲线方程为：该解也为一精确解。3.4四边简支矩形薄板受均布荷载q，求挠曲方程。微分方程：D4wq=0，yyxBOACab图3.6四边简支矩形薄板边界条件：。取满足边界条件的一级近似试函数：残值：(1)连续型方法由甲醛积分式解得：对于方板（a=b），中心挠度wmax=0.004146qa4/D，精确解wmax=0.004168qa4/D，误差2.5%。(2)离散型方法取均匀内部配点，将薄板分成NN网格。若N=4，则有9个配点，如图3.7所示。AAyxBOC图3.7板内配点取法求出各点残值Ri（i=1,2,…,9）（Ri是的线性表达式），并写成矩阵形式{R}=[C]nm{}{b}残值平方和：Id={R}T{R}={}T[C]T[C]{}{b}T[C]{}{}T[C]T{b}+{b}T{b}由（标量对向量的导数为零），得线性代数方程：[C]T[C]{}=[C]T{b}。本例中：，代入上式，解得解得：当网格数目增大，值收敛于连续型的结果。(3)四边简支矩形薄板中一点(,)处作用一竖向集中荷载P，求板的挠曲方程和中心挠度。yyxBOACab图3.8四边简支矩形薄板受集中荷载作用微分方程：D4wP(x)(y)=0，取满足边界条件的一阶近似试函数：残值：加权残数方程：由函数的性质：，可解得习题根据试函数分类，加权残数法可分为哪几种方法？写出相应的加权积分式的一般表达式，并说明试函数的一般选取原则。根据权函数分类（采用内部法），加权残数法可分为哪几种基本方法？其权函数和加权积分式各是什么？四边简支矩形薄板，边长分别为a和b，弯曲刚度为D。试用伽辽金方法求板在均布荷载q作用下的挠曲方程和中心挠度（取一阶近似）。矩形薄板边长为a和b，弯曲刚度为D，四边固支。试用伽辽金方法求板在均布荷载q作用下的挠曲方程和中心挠度（取一阶近似）。四边简支的矩形薄板，边长分别为a和b。试用连续型最小二乘法求板在均布荷载q作用下的挠曲方程和中心挠度（一阶近似）。图示四边简支矩形薄板，边长分别为a和b。板在x=，y=处作用一竖向集中荷载P，试用连续型伽辽金方法求板的挠曲方程和中心挠度（取一阶近似）。图示四边简支矩形薄板，边长分别为a和b。已知板在x=线上作用均布线荷载p，试用连续型最小二乘法求板的挠曲方程和中心挠度（取一阶近似）。yxyxzpabyxzPab题3.6图题3.7图

第4章边界单元法4.1主要特点及发展概况4.1.1边界单元法与有限单元法的优点：应用范围广，功能强大；但也有其弱点：(1)如单元总数一般要求较多，计算工作量大，尤其对三维问题；(2)单元之间的有些物理量可能不连续（如用位移法，位移连续，但应力应变不一定连续），出现不合理的跳跃间断；(3)不易处理（半）无限区域问题、应力集中问题。边界单元法的主要优点：(1)只需将结构的边界进行离散，降低了问题的维数，使计算工作量大大降低，计算准备工作也大为简化；(2)计算精度相对较高，误差仅来源于边界；(3)能较好处理应力集中问题、无限、半无限区域问题；(4)可与有限元法、差分法等其他数值方法联合应用，解决复杂的工程技术问题。边界单元法的不足：需要较多的数学知识，不如有限单元法那么简单直观，易于接受；最后得到的代数方程组的系数矩阵是个非对称满阵；不如有限单元法成熟以及具有广泛的应用范围。4.1.2边界单元法的1903年Fredholm提出了将由微分方程描述的边值问题变换为边界积分方程描述，并完成了势问题基本方程的变换，证明了解答的存在性和唯一性，但无法获得其解析解。1963年Jaswon和Symm提出了求解Fredholm边界积分方程的数值解法；1965年Kupradze对弹性力学边值问题列出了向量积分方程；1968~69年Cruse给出了三维弹性力学积分方程的数值解；1972年Cruse这种边界积分方程方法应用到断裂力学、塑性力学问题中；1978年Brebbia在其著作“TheBoundaryElementMethodforEngineering”中首次引入了“边界单元法”的术语，并阐述了其基本原理及应用。其后边界单元法的研究及应用得到了更迅速的发展。4.2基本原理——二维势问题的边界单元法基本思想：首先将问题的微分方程变换为边界积分方程，再对边界进行离散，形成关于边界离散物理量的线性代数方程组，解出边界值，域内值最后通过边界积分求出。4.2.1边界积分方程的形成——直接法公式变换方法：Green公式、加权残数法、虚功原理。势问题的基本方程可用Laplace方程描述：，其中xx1ox212图4.1二维势问题示例(1)域内点的边界积分方程应用加权残数法，选择满足边界条件的试函数u，得到加权残数方程，其中W为权函数。利用Green公式，令=W，=u，则有（通过两次分步积分也可获得）：取权函数满足下列方程式，满足该方程的解称为基本解。如对温度场问题，该基本解的物理意义就是无限大物体上一点处有一单位热源作用时，引起的在场内x点上的温度。W表示源点与场点x间两点状态的影响函数。对上述二维Laplace问题，基本解为根据函数的性质：有于是可得：(a)该式将域内点的对应值u()用边界积分方程表示出来。(2)边界点的边界积分方程若将上式直接用于边界点的计算，则因在x处（r0）出现奇异积分，故其中的积分可能出现跳跃。为此，将奇异积分分离出来计算，即=(-)+，n图4.2奇异积分区域的分离当半径0时，考虑到，故在上有：于是当0时，-=，因此光滑边界上任意点处的值可用边界积分方程表示为(3)边界单元法精确求解上述边界积分比较困难，因此采用对边界进行离散的数值方法——边界单元法求解。求解步骤：二