现代接触动力学2.doc

现 代 接 触 力 学
相应的偏微分方程必须在整个区域上得到满足， 见图3．5

物体的边界分为两类和，在边界上位移受到约束，在边界上应力受到约束。在上的边界条件称之为位移边界条件，或强制边界条件，或Neumann边界条件。对大多数实用问题，位移边界条件通常为
（3．47）
在上的边界条件称之为力的边界条件，或自然边界条件，或Neurmann边界条件。力的边界条件可写为
（3．48）
现以偏微分方程组式（3．43）为例，说明近似解的构造。位移由个待定参数和一族已知函数来近似，即
（3．49）
其中称之为插值函数，或形函数，或基函数，或试探函数；为待定参数，与结点在方向的位移有关。用矩阵形式，式（3．49）可表达为
（3．50）
其中
（3．51）
（3．52）
（3．53）
形函数的合理选择将在第（3．6）节中详细讨论。由应力平衡方程式（3．41）可获得余量
（3．54）
等效积分的加权余量形式为
（3．55）
通过合理地选择个线性独立的权函数，式（3．50）中的个待定参数可由式（3．55）唯一确定。对于精确解，余量恒等于零。而对于近似解，余量则不总是为零，式（3．55）表明余量在权函数方向投影的积分为零。将余量表达式（3．54）代入式（3．55）中，则有
（3．56）
通过乘积函数的微分法则
（3．57）
和分部积分，式（3．56）可转化为
（3．58）
利用高斯积分定理（见附录式（A．6．8）），得
（3．59）
其中为域的边界，为边界的外法线方向。将式（3．59）代入式（3．58），则体积分
（3．60）
在位移边界上的积分为零，因为位移边界条件应被满足，权函数在位移边界上选为零，如果将权函数看成容许的变分，则在位移边界上变分必须为零，即。在力的边界上
（3．61）
代入式（3．60），则有
（3．62）
推导式（3．62）的另一种途径是对力的边界条件式（3．48）同样引入权函数，要求权函数对偏微分方程余量的区域积分和对力的边界条件余量的边界积分之和为零。这时虽然有更多的积分项，但部分积分项正好相互抵消，见文献[Papadopoulos 1996a].。
在等效积分形式（3．56）中，位移必须满足位移边界条件和力的边界条件。而在等效积分的弱解形式（3．62）中，力的边界条件已被考虑，位移只需满足位移边界条件。鉴于使用整体积分和加权余量，力的平衡方程和边界条件只需近似地满足。由于在分部积分中通过提高权函数的可微性要求而降低了对位移函数的二阶可微性要求，积分形式因此被称之为弱解形式或变分形式。具体地说，位移函数在弱解形式中只要一次可微即可，不再需要二次可微。由于权函数通常相对简单一些 ，通过合理的选择，权函数的可微性要求容易满足。
为了获得确定待定系数的线性方程，将近似的位移场式（3．50）代入应变表达式（3．19）中，则得应变张量
（3．63）
相应的应变和位移间的矩阵表达式为
（3．64）
其中
（3．65）
（3．66）
将式代入材料的本构关系式中，则有
（3．68）
（3．69）
（3．70）
（3．71）
以及引进与矩阵和结构类似的矩阵
（3．72）
（3．73）
弱解形式（3．62）可写成
（3．74）
由本够关系式（3．67），则进一步有
（3．75）
由于待定参数与空间坐标无关，可以提到积分号外，因此式（3．75）可写为
（3．76）
于是，所求的线性方程组为
（3．77）
其中
（刚度矩阵） （3．78）
（质量矩阵） （3．79）
（体力向量） （3．80）
（面力向量） （3．81）
方程式（3．77）是一组线性常微分方程，给定初始条件，可用第7章7．2节叙述的数值积分方法求解。对于静力学问题，因此方程式（3．77）为一组线性代数方程
（3．82）
求解线性代数方程组的基本算法参见第7章7．5节。不论是静力学问题还是动力学问题，求解出待定参数后，位移项即可由式（3．50）确定。
方程式（3．77）涉及整个区域，为使微分方程的近似解有满意的精度，人们必须适当地选取权函数和形函数。对复杂的几何形体，它们的确定通常是一项困难的工作。在有限元法中，区域有许多简单的微小单元，即有限元来近似。其优点是对简单的有限单元权函数和形函数很容易确定，然后通过合理的装配，偏微分方程在整个区域上的解则能够近似地计算。权函数的不同选择导致不同的计算方法，本章3．3节简要地讨论了权函数的若干选择。权函数的选择取决于插值的精度、单元的几何形状以及待定系数的物理意义等因素，这些将在本章3．6节中讨论。待定参数的物理意义取决于应用的单元类型。对线弹性问题，用有限元法求解位移场时，这些待定参数为结点的位移。值得注意的是，到目前为止边界条件中只涉及力的边界条件；位移边界条件只有在求解方程式（3．77）中才被考虑（参见本章节）
本节中方程式由加权余量法导出，更详细的描述可参看[Burnett 1987]。人们亦可利用变分法，例如瑞利-迦辽金法来推导有限元方程，如在参考书[Reddy 1993]中那样。这种方法由于类似于结构力学中的极值原理，通常受到人们的偏爱，但在应用中变分法对偏微分方程的要求甚严，且不易被理解和应用。
3．3 加权余量法
为使方程组式（3．82）中的待定系数能够唯一确定，权函数必须满足某些要求，诸如线性独立性。在文献中，人们引进了各种不同的权函数，参见[Finlayson Scriven 1966a]。本书采用迦辽金法，取权函数为形函数，即
（3．83）
其优点是质量矩阵和刚度矩阵K为对称矩阵。在某些情况下，迦辽金法和瑞利-里兹法导致同样的结果，参见[Reddy 1993]。
配点法的基本思想是微分方程只要求在若干点上得到满足。权函数相当于取狄拉克函数。该方法强迫余量在区域个点上等于零。配点法的精度与所选的配点位置密切相关。一般来说配点法中的刚度矩阵和质量矩阵为非对称矩阵。配点法的一种改进方法是子域法，这种方法强迫余量在子域上的积分为零。权函数在子域内为1，在子域外为零。
最小二乘法的实质是使得余量平方和取最小值。它可由变分原理而导出，参见[Papadopoulos 1996a]。其权函数。最小二乘法中的刚度矩阵和质量矩阵同样为对称矩阵。
3．4 基本边界条件
对位移边界条件的处理在[Rurnett 1987]一书中有着详细的描述，这里仅讨论有关算法和程序实现的若干细节。位移边界条件可以采用以下方法引入。
1） 矩阵分块法
在用有限元法求解位移场中，位移边界上位移已知，力的边界上力已知，这样在方程
（3．84）
中，每一行要么结点位移已知，要么结点力已知，为了求解未知的变量，方程式（3．84）需要进行适当变动。按结点，位移是已知还是待定需要重新组合方程。
通过调换矩阵和列向量的行以及矩阵的行和列，方程式（3．84）可用分块矩阵形式表示，即
（3．85）
所有在位移边界上给定的位移变量包含在中，未知的位移变量包含在中，所有在力的边界上给定的力包含在中，未知的反力包含在中。利用分块矩阵方程的第一块，由下式
（3．86）
（3．87）
（3．88）
可求得未知的结点位移。
在数值求解中，矩阵的逆无需直接计算，而是采用在第7章7．5节描述的方法求解线性方程组式（3．87）。在解出结点位移后，位移边界上的约束反力可由简单的乘法运算而确定，即
（3．89）
在许多有限元程序中刚度矩阵并未建立，实际使用的只是子矩阵，作用力则通过和作用应的修正。子矩阵和只是在需要计算约束力时才建立。矩阵分块法虽然原理简单，但其数据结构的操作却相当费时，因为结点位移的顺序性被破坏。
2） 修改元素法
在修改元素法中，结点位移的顺序无需改变，但矩阵和力向量中的元素须作相应的修改。这里简要地叙述其操作过程：对方程组（3．84）进行处理。首先检验在第行中是否结点位移已给定，如是，则对力向量中的每个元素分别去掉力的修正项。在符号上所有未知的结点力则写成。对每一行进行上述操作后，则形成非对称矩阵，修正的力向量和未知变向量。线性方程组
（3．90）
则可直接求解。方程组的解既包括未知的解点位移又包括未知的结点约束反力。
现借助一个简单的例子说明这一方法。原线性方程组为
（3．91）

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

图3.6 位移边界条件

如果平面如图3．6（b）所示，倾斜角为，则自由和约束的位移方向不再平行于坐标轴。引进新的坐标系，则待求的量为和，已知的量为和。为表达这种坐标变换，方程式（3．84）中与点有关的行和列须乘以变换矩阵
（3．92）
其他所有的行和列则保持不变。借助于正交变换矩阵
（3．93）
方程式（3．84）在形式上可转换为
（3．94）
由于矩阵的稀疏性，在上述表达式的数值计算中，自然用不着进行矩阵乘法运算，而只需对有关的行和列作相应的修正。
在坐标系中算出的量及通过变换式及即可转换为在坐标系中相应的量。通过合理的选择变换矩阵，人们能够方便地处理位移变量间的线性相关性，即
（3．95）
一般来说变换矩阵
（3．96）
不再具有对称性和正交性。比如，人们要求结点和具有相同的位移，则可转换为
（3．97）
在新的坐标系中，则要求位移边界条件必须得到满足，得到解后，通过变换式则有
（3．98）
（3．99）
可见，结点和总有相同的位移。对有些边界条件，比如，人们可以称对地处理，例如要求。然而引进这些约束方程时，人们必须注意约束方程间的独立性和相容性。
3．6 单元格式
根据式（3．50），位移场由插值函数近似描述。由此产生的误差则取决于单元的选择以及插值函数的类型和阶次。下面将对不同类型的单元作一简要描述。
3．6．1 等参单元
依据文献[Brenner Scott 1994]，一个有限单元由三个要素来定义。如果：
（1） 为一个分段光滑的有界区域（单元区域）；
（2） 为在域上的一有限维函数空间（插值函数空间）；
（3） 为函数空间的基（结点变量）。
则称为一有限单元。
在有限元中经常使用等参单元。所谓等参单元是指单元几何形状
（3．100）
和单元的位移场函数
（3．101）
采用相同的插值函数来近似描述，其选择虽然十分自然，但并非必须如此。在非等参单元中，人们可以对单元几何形状和位移场分别采用不同的插值函数来近似描述。
3．6．2 分类

依据文献[Dhatt Touzot 1983]，有限单元有下述不同的特征：
（1） 形状（比如三角形单元，四边形单元）；
（2） 个插值点的自然坐标；
（3） 自由度个数；
（4） 结点坐标的定义（拉格朗日单元：位移，Hermite单元：位移和斜率）；
（5） 多项式插值函数的基（比如线性，二次，完整，非完整）；
（6） 单元在结点处的连续性（比如，）符号表示所有至阶连续可微的函数空间。
表3．1列出了5种平面单元，它们皆为拉格朗日单元，在结点具有连续性。
表31 若干单元描述
单元 结点数目 自由度 多项式 基
三角形 3 6 线性
三角形 6 12 二次
四边形 4 8 双线性
四边形 8 16 二次
四边形 9 18 二次

3．6．3 线性三角形单元
对大多数单元来说，一种行之有效的方法是建立参考单元，然后将真实单元和参考单元通过一对一对的映射联系起来。这种方法不仅直观而且还可以避免不必要的重复运算。图3．7演示在参考区域平面内三角形单元和它在真实区域平面内的映射。参考单元由不等式确定。通过映射
（3．102）
可确定参考单元相对应的真实单元。
人们亦可把坐标系看作为描述真是单元的非直角局部自然坐标系。对三角形单元，由式（3．100）和（3．101）可得
（3．103）
其中线性插值函数需满足下述条件：
（3．104）
表3．2列出了这些插值函数及它们对自然坐标的导数。
表3．2 线性三角形单元插值函数及其导数

三角形单元的插值函数被描绘在图3．8中，他们具有下述重要特征：
（1） 参考单元的结点被映射到真实的结点

（2） 参考单元的边界被映射到真实单元的边界。例如由确定的边界线段
段在真实区域内为
（3．105）
线段经过线性映射后仍为一条线段，所有的点位于结点和之间。
（3） 映射必须一一对应，因此雅克比矩阵 （3．106）
在整个单元区域内必须非奇异。在实际计算中，雅克比矩阵由下式计算：
（3．107）
相应的元素为
（3．108）
（3．109）
（3．110）
（3．111）
对线性三角形单元有
（3．112）
如果三个结点不位于一条直线上，则行列式。由此矩阵非奇异，映射为一一对应。如果三个结点按逆时针排序，则行列式。
应变矩阵可由导数按式（3．66）计算，而导数则由表3．2中给出的插值函数对自然坐标的导数及雅克比矩阵共同确定，即
（3．113）
因为雅克比矩阵非奇异，，所以逆矩阵在整个单元内存在，人们可以确定应变矩阵以及计算刚度矩阵（3．78）中的求积项。
刚度矩阵易可在参考单元区域内求积而得，若用表示真实单元区域，表示参考单元区域，则有关系式
（3．114）
类似地对质量矩阵有关系式
（3．115）
值得注意的是，插值函数及其导数与真实单元的结点坐标无关，因此一些常量