现代接触动力学2.doc

1、现 代 接 触 力 学相应的偏微分方程必须在整个区域上得到满足， 见图 35o 1x2x物体的边界分为两类 和 ，在边界 上位移受到约束，在边界 上应力受到约束。uquq在 上的边界条件称之为位移边界条件，或强制边界条件，或 Neumann 边界条件。对大多u数实用问题，位移边界条件通常为()uuXX常 数（347）在 上的边界条件称之为力的边界条件，或自然边界条件，或 Neurmann 边界条件。力q的边界条件可写为()qTXnXA（348）现以偏微分方程组式（343）为例，说明近似解的构造。位移 由 个待定参数iun和一族已知函数 来近似，即()iU()jN(1,)ijijuXUn（349

2、）其中 称之为插值函数，或形函数，或基函数，或试探函数； 为待定参数，与()jX ()ijU结点在 方向的位移有关。用矩阵形式，式（349）可表达为iuNUA（350）其中 12uR（351）12 200nnNN （352）现 代 接 触 力 学 211212212()()()()nnUUUR （353）形函数的合理选择将在第（36）节中详细讨论。由应力平衡方程式（341）可获得余量()(,12)iki iTRubukx（354）等效积分的加权余量形式为()0(1,2;,)judijn w（355）通过合理地选择 个线性独立的权函数 ，式（350）中的 个待定参数可由式njw2n（355）唯一

3、确定。对于精确解，余量 恒等于零。而对于近似解，余量则不总是为()iRu零，式（355）表明余量在权函数方向投影的积分为零。将余量表达式（354）代入式（355）中，则有()0ikjiiTbudxw（356）通过乘积函数的微分法则()()()aabbxbx（357）和分部积分，式（356）可转化为()()0j jikik jiiTTddbudxxww（358）利用高斯积分定理（见附录式（A68） ） ，得()ijikjkjiqdTndqx（359）其中 为域 的边界， 为边界的外法线方向。将式（359）代入式（358） ，则体积分()0jikjijijiTdqbudxWW（360）在位移边界上

4、的积分 为零，因为位移边界条件应被满足，权函数在位移边界uji现 代 接 触 力 学上选为零，如果将权函数 看成容许的变分 ，则在位移边界上变分必须为零，即ujWiu。在力的边界 上()0,i uXqnTA（361）代入式（360） ，则有()qjikjijijiTdbudxWW（362）推导式（362）的另一种途径是对力的边界条件式（348）同样引入权函数，要求权函数对偏微分方程余量的区域积分和对力的边界条件余量的边界积分之和为零。这时虽然有更多的积分项，但部分积分项正好相互抵消，见文献Papadopoulos 1996a.。在等效积分形式（356）中，位移必须满足位移边界条件和力的边界条件

5、。而在等效积分的弱解形式（362）中，力的边界条件已被考虑，位移只需满足位移边界条件。鉴于使用整体积分和加权余量，力的平衡方程和边界条件只需近似地满足。由于在分部积分中通过提高权函数的可微性要求而降低了对位移函数的二阶可微性要求，积分形式因此被称之为弱解形式或变分形式。具体地说，位移函数 在弱解形式中只要一次可微即可，不再需要u二次可微。由于权函数 通常相对简单一些 ，通过合理的选择，权函数的可微性要求容jW易满足。为了获得确定待定系数的线性方程，将近似的位移场式（350）代入应变表达式（319）中，则得应变张量()()11()22()j jki ikijjijijk jkik ikji ji

6、uNUxxxNU（363）相应的应变和位移间的矩阵表达式为BUA（364）其中现 代 接 触 力 学3112R（365）1211 3222211100nnnNNXXBR （366）将式代入材料的本构关系式中，则有121 1122 11 0jJJJj JJJTTxXxxAWW（368）11220jjjbbAW（369）11220jjjqqA（370）112200jj jj juuNUAAWW（371）以及引进与矩阵 和 结构类似的矩阵BN现 代 接 触 力 学1211 3222211100nnnxxxRxxx WW（372）12 200nnRWW（373）弱解形式（362）可写成 qTTTTrd

7、NbdNUddAAAAWWWW（374）由本够关系式（367） ，则进一步有 qTTTTrCUdNbdNddAAAAWWWW（375）由于待定参数与 空间坐标无关，可以提到积分号外，因此式（375）可写为 qTTTTrBCdNbddNdUAAAAWWW（376）于是，所求的线性方程组为rKUfMAA（377）其中（刚度矩阵） 2TnKBCdRAW（378）（质量矩阵） 2nMN现 代 接 触 力 学（379）（体力向量） 2TnfNbdRAW（380）（面力向量） 2qTnrf（381）方程式（377）是一组线性常微分方程，给定初始条件，可用第 7 章 72 节叙述的数值积分方法求解。对于静力

8、学问题 ，因此方程式（377）为一组线性代数方程0urKUfA（382）求解线性代数方程组的基本算法参见第 7 章 75 节。不论是静力学问题还是动力学问题，求解出待定参数 后，位移项即可由式（350）确定。u方程式（377）涉及整个区域 ，为使微分方程的近似解有满意的精度，人们必须适当地选取权函数 和形函数 。对复杂的几何形体 ，它们的确定通常是一项困难的NW工作。在有限元法中，区域 有许多简单的微小单元，即有限元来近似。其优点是对简单的有限单元权函数 和形函数 很容易确定，然后通过合理的装配，偏微分方程在整个区域 上的解则能够近似地计算。权函数 的不同选择导致不同的计算方法，本章NW33

9、节简要地讨论了权函数的若干选择。权函数的选择取决于插值的精度、单元的几何形状以及待定系数的物理意义等因素，这些将在本章 36 节中讨论。待定参数的物理意义取决于应用的单元类型。对线弹性问题，用有限元法求解位移场时，这些待定参数为结点的位移。值得注意的是，到目前为止边界条件中只涉及力的边界条件 ；位移边界0TnqA条件 只有在求解方程式（377）中才被考虑（参见本章节）0u本节中方程式由加权余量法导出，更详细的描述可参看Burnett 1987。人们亦可利用变分法，例如瑞利-迦辽金法来推导有限元方程，如在参考书Reddy 1993中那样。这种方法由于类似于结构力学中的极值原理，通常受到人们的偏爱

10、，但在应用中变分法对偏微分方程的要求甚严，且不易被理解和应用。33 加权余量法为使方程组式（382）中的待定系数 能够唯一确定，权函数 必须满足某些要求，ujW现 代 接 触 力 学诸如线性独立性。在文献中，人们引进了各种不同的权函数，参见Finlayson Scriven 1966a。本书采用迦辽金法，取权函数为形函数，即,jjNBW（383）其优点是质量矩阵 和刚度矩阵 K 为对称矩阵。在某些情况下，迦辽金法和瑞利-里兹法M导致同样的结果，参见Reddy 1993。配点法的基本思想是微分方程只要求在若干点上得到满足。权函数 相当于取狄拉克j函数 。该方法强迫余量在区域 个点 上等于零。配点

11、法的(),()jjjxx即 Wnjx精度与所选的配点位置 密切相关。一般来说配点法中的刚度矩阵 和质量矩阵 为非j KM对称矩阵。配点法的一种改进方法是子域法，这种方法强迫余量在 子域上 的积分为j零。权函数 在子域 内为 1，在子域外 为零。jjj最小二乘法的实质是使得余量平方和取最小值。它可由变分原理而导出，参见Papadopoulos 1996a。其权函数 。最小二乘法中的刚度矩阵 和质量矩阵jjRuWK同样为对称矩阵。M34 基本边界条件对位移边界条件的处理在Rurnett 1987一书中有着详细的描述，这里仅讨论有关算法和程序实现的若干细节。位移边界条件可以采用以下方法引入。1） 矩

12、阵分块法在用有限元法求解位移场中，位移边界 上位移已知，力的边界 上力已知，这样uq在方程KUfA（384）中，每一行要么结点位移 已知，要么结点力 已知，为了求解未知的变量，方程式iuif（384）需要进行适当变动。按结点，位移是已知还是待定需要重新组合方程。现 代 接 触 力 学通过调换矩阵 和列向量 的行以及矩阵 的行和列，方程式（384）可用分块矩UfK阵形式表示，即ququKfA（385）所有在位移边界上给定的位移变量包含在中，未知的位移变量包含在中，所有在力的边界上给定的力包含在中，未知的反力包含在中。利用分块矩阵方程的第一块，由下式qquqKUfA（386）qqufAA（387）

13、1()qquUKfUA（388）可求得未知的结点位移 。q在数值求解中，矩阵 的逆无需直接计算，而是采用在第 7 章 75 节描述的方法求qK解线性方程组式（387） 。在解出结点位移 后，位移边界上的约束反力 可由简单的qUuf乘法运算而确定，即11()uquuqququufKUKfKUAAAA（389）在许多有限元程序中刚度矩阵 并未建立，实际使用的只是子矩阵 ，作用力 则qKqf通过 和 作用应的修正。子矩阵 和 只是在需要计算约束力时才建立。矩阵分quKuqK块法虽然原理简单，但其数据结构的操作却相当费时，因为结点位移的顺序性被破坏。2） 修改元素法现 代 接 触 力 学在修改元素法中

14、，结点位移的顺序无需改变，但矩阵 和力向量 中的元素须作相应Kf的修改。这里简要地叙述其操作过程：对方程组（384）进行处理。首先检验在第 行中i是否结点位移 已给定，如是，则对力向量中的每个元素 分别去掉力的修正项 。iUif ijKU在符号上所有未知的结点力 则写成 。对每一行进行上述操作后，则形成非对称矩阵ifiU，修正的力向量 和未知变向量 。线性方程组KfKfA（390）则可直接求解。方程组的解 既包括未知的解点位移又包括未知的结点约束反力。U现借助一个简单的例子说明这一方法。原线性方程组为112233445092KfUfA（391） 1234509124fUfA1234105.15

15、399.24UfKff A运用修改元素法，其解则为 1234,10fuf修改元素法尽管看上去完美，但由于矩阵 的非对称性，在计算中并不比其他方法有K效，因而很少使用。现 代 接 触 力 学在上述两种引入位移边界条件的方法中，方程组需要重新组合，然而可以获得在位移边界上准确的约束反力。人们亦可采用其他途径，由插值函数的导数和结点位移获得约束反力的近似解，参见Reddy 1993。这种方法不用重新组合方程组，因此能够节省计算时间，但以近似误差为代价。如果在时间积分中需要约束反力，则最好运用准确计算约束反力的方法。35 坐标变换通常坐标轴的方向选取与结点位移平行的方向。对图 36（a）中所示的在水平

16、面上滑动的结点 ，其边界条件为： 待求， 已知， 待求。Pxuxfyf()b()aPP图 3.6 位移边界条件如果平面如图 36（b）所示，倾斜角为 ，则自由和约束的位移方向不再平行于坐标轴。引进新的坐标系，则待求的量为 和 ，已知的量为 和 。为表达这种坐标xuyf yuxf变换，方程式（384）中与 点有关的行和列须乘以变换矩阵Pcosini（392）其他所有的行和列则保持不变。借助于正交变换矩阵1cosini1P （393）现 代 接 触 力 学方程式（384）在形式上可转换为,TTTKPUfPfKUfAAA（394）由于矩阵 的稀疏性，在上述表达式的数值计算中，自然用不着进行矩阵乘法运

17、算，而只需对有关的行和列作相应的修正。在坐标系 中算出的量 及 通过变换式 及 即可转换为在,xyfAfP坐标系 中相应的量。通过合理的选择变换矩阵 ，人们能够方便地处理位移变量间, P的线性相关性，即ij iaUaUgk（395）一般来说变换矩阵1 j11ji i iaaiPa 第 行第 行第 行kk（396）不再具有对称性和正交性。比如，人们要求结点和具有相同的位移，则可转换为0ijiU（397）在新的坐标系中，则要求位移边界条件 必须得到满足，得到解 后，通过变换i jU式 则有UPA1jjjjPUA（398）0iijijijjpU（399）可见，结点和总有相同的位移。对有些边界条件，比

18、如， 人们可以称对地,ijk现 代 接 触 力 学处理，例如要求 。然而引进这些约束方程时，人们00ijiiUU和 k必须注意约束方程间的独立性和相容性。36 单元格式根据式（350） ，位移场 由插值函数近似描述。由此产生的误差则取决于单元的选择以及插值函数的类型和阶次。下面将对不同类型的单元作一简要描述。361 等参单元依据文献Brenner Scott 1994，一个有限单元由三个要素来定义。如果：（1） 为一个分段光滑的有界区域（单元区域） ；KR（2） 为在域 上的一有限维函数空间（插值函数空间） ；P（3） 为函数空间 的基（结点变量） 。12,N kP则称 为一有限单元。在有限元

19、中经常使用等参单元。所谓等参单元是指单元几何形状xNA（3100）和单元的位移场函数uU（3101）采用相同的插值函数来近似描述，其选择虽然十分自然，但并非必须如此。在非等参单元中，人们可以对单元几何形状和位移场分别采用不同的插值函数来近似描述。362 分类依据文献Dhatt Touzot 1983，有限单元有下述不同的特征：（1） 形状（比如三角形单元，四边形单元） ；（2） 个插值点的自然坐标 ；nn（3） 自由度个数； （4） 结点坐标的定义（拉格朗日单元：位移，Hermite 单元：位移和斜率） ；（5） 多项式插值函数的基（比如线性，二次，完整，非完整） ；（6） 单元在结点处的连续

20、性（比如 ， ）符号 表示所有至 阶连续可微的01,CsS函数空间。表 31 列出了 5 种平面单元，它们皆为拉格朗日单元，在结点具有 连续性。0C表 31 若干单元描述 单元 结点数目 自由度 多项式 基三角形 3 6 线性 1,现 代 接 触 力 学三角形 6 12 二次 21,四边形 4 8 双线性 |四边形 8 16 二次 22,四边形 9 18 二次 21363 线性三角形单元对大多数单元来说，一种行之有效的方法是建立参考单元，然后将真实单元和参考单元通过一对一对的映射联系起来。这种方法不仅直观而且还可以避免不必要的重复运算。图37 演示在参考区域 平面内三角形单元和它在真实区域平面

21、内的映射。参考单元 由 rV不等式 确定。通过映射01,123():();,)reVxx（3102）可确定参考单元相对应的真实单元 。eV人们亦可把坐标系 看作为描述真是单元 的非直角局部自然坐标系。对三角形,e单元，由式（3100）和（3101）可得112323() ,xxNINIuNU AAA（3103）其中线性插值函数 需满足下述条件：10ijijN当当（3104）表 32 列出了这些插值函数及它们对自然坐标的导数。表 32 线性三角形单元插值函数及其导数/1112 030i i iiNN现 代 接 触 力 学三角形单元的插值函数被描绘在图 38 中，他们具有下述重要特征：（1） 参考单

22、元的结点被映射到真实的结点123100xxIA结 点 ： ，123xxI结 点 2： ，123010xxIA结 点 3： ，（2） 参考单元的边界被映射到真实单元的边界。例如由 确定的边界线段,0段 在真实区域内为11221310xxIxA（3105）线段 经过线性映射后仍为一条线段，所有的点位于结点 和 之间。12 1x2（3） 映射必须一一对应，因此雅克比矩阵 Jy（3106）在整个单元区域内必须非奇异。在实际计算中，雅克比矩阵 由下式计算：JiiNxJ（3107）相应的元素为1iixJ（3108）现 代 接 触 力 学12iiNxJ（3109）21iiyJ（3110）2iiNyJ（311

23、1）对线性三角形单元有213121313121,detxxJJxyxyyy（3112）如果三个结点不位于一条直线上，则行列式 。由此矩阵 非奇异，映射为一一对et0JJ应。如果三个结点按逆时针排序，则行列式 。应变矩阵 可由导数 按式（366）计算，而导数 则由表 32 中给B/iNx/iNx出的插值函数对自然坐标的导数 及雅克比矩阵 共同确定，即J1 1,iii iiijjjJxx AA即 kk（3113）因为雅克比矩阵 非奇异， ，所以逆矩阵 在整个单元内存在，人们可以确定Jdet0J1J应变矩阵 以及计算刚度矩阵（378）中的求积项 。BTBCA刚度矩阵 易可在参考单元区域内求积而得，若用 表示真实单元区域， 表示参Ker考单元区域，则有关系式dte rTeTrBCdJA（3114）类似地对质量矩阵 有关系式Mete rTeTrNNJA（3115）值得注意的是，插值函数 及其导数 与真实单元的结点坐标无关，因此一些常量i /i