现代接触动力学2

相应的偏微分方程必须在整个区域上得到满足， 见图 35 o 1 x 2 x 物体的边界分为两类和，在边界上位移受到约束，在边界上应力受到约束。 u q u q 在上的边界条件称之为位移边界条件，或强制边界条件，或 Neumann 边界条件。对大多 u 数实用问题，位移边界条件通常为 () u u XuX常数 （347） 在上的边界条件称之为力的边界条件，或自然边界条件，或 Neurmann 边界条件。力 q 的边界条件可写为 ()() q T Xn XqXA （348） 现以偏微分方程组式（343）为例，说明近似解的构造。位移由个待定参数 i un 和一族已知函数来近似，即() i U() j NX ()()(1, ) ijij uNX Ujn （349） 其中称之为插值函数，或形函数，或基函数，或试探函数；为待定参数，与() j NX() ij U 结点在方向的位移有关。用矩阵形式，式（349）可表达为j i X uN UA （350） 其中 12 1 2 u uR u （351） 122 2 12 000 000 nn n NNN NR NNN （352） 21 1 12 11 22212 ()()()()()() n nn UUUUUUUR （353） 形函数的合理选择将在第（36）节中详细讨论。由应力平衡方程式（341）可获得余量 ( )( ,1,2) ik ii k T R ubui k x （354） 等效积分的加权余量形式为 ( )0(1,2;1, ) jR u d ijn w （355） 通过合理地选择个线性独立的权函数，式（350）中的个待定参数可由式n j w2n （355）唯一确定。对于精确解，余量恒等于零。而对于近似解，余量则不总是为( ) i R u 零，式（355）表明余量在权函数方向投影的积分为零。将余量表达式（354）代入式 （355）中，则有 ()0 ik jii k T bu d x w （356） 通过乘积函数的微分法则 ( )( ( ) ( )( ) ( )( ) a xa x b xb x b xa x xxx （357） 和分部积分，式（356）可转化为 () ()0 jjik ikjii kk T T ddbu d xx ww w （358） 利用高斯积分定理（见附录式（A68） ） ，得 () i jik jikkji k q T dT n dq d x w ww （359） 其中为域的边界，为边界的外法线方向。将式（359）代入式（358） ，则体积n 分 ()0 j ikjijiji k T dq dbu d x W WWW （360） 在位移边界上的积分为零，因为位移边界条件应被满足，权函数在位移边界 u ji q d W 现 代 接 触 力 学 上选为零，如果将权函数看成容许的变分，则在位移边界上变分必须为零，即 u j W i u 。在力的边界上()0, iu u XX q qn TqA （361） 代入式（360） ，则有 () q j ikjijiji k T dbu dq d x W WWW （362） 推导式（362）的另一种途径是对力的边界条件式（348）同样引入权函数，要求权函 数对偏微分方程余量的区域积分和对力的边界条件余量的边界积分之和为零。这时虽然有 更多的积分项，但部分积分项正好相互抵消，见文献Papadopoulos 1996a.。 在等效积分形式（356）中，位移必须满足位移边界条件和力的边界条件。而在等效积分 的弱解形式（362）中，力的边界条件已被考虑，位移只需满足位移边界条件。鉴于使用 整体积分和加权余量，力的平衡方程和边界条件只需近似地满足。由于在分部积分中通过 提高权函数的可微性要求而降低了对位移函数的二阶可微性要求，积分形式因此被称之为 弱解形式或变分形式。具体地说，位移函数在弱解形式中只要一次可微即可，不再需要u 二次可微。由于权函数通常相对简单一些 ，通过合理的选择，权函数的可微性要求容 j W 易满足。 为了获得确定待定系数的线性方程，将近似的位移场式（350）代入应变表达式 （319）中，则得应变张量 () ()11 ()() 22 11 ()() 22 jjkk iikk ij jiji jkjk ikik kkk jiji uN U uN U xxxx NN NN UUU xxxx （363） 相应的应变和位移间的矩阵表达式为 B UA （364 ） 其中 3 1 112212 2R （365） 12 111 3 2 12 222 1122 212121 000 000 n n n nn NNN XXX NNN BR XXX NNNNNN XXXXXX （366） 将式代入材料的本构关系式中，则有 112111 1212 22 121222 21 12 0 0 j JJJ jJJ J TTT xXxx T TTT xx xX A W WWW WWW W （368） 1 1 2 2 0 0 jj jj b bb bb A WW WW （369） 1 1 22 0 0 j j j j q q q qq A WW WW （370） 1 1 22 00 00 j jj jj j u u N U uu AAA WWW WWW （371） 以及引进与矩阵和结构类似的矩阵BN 现 代 接 触 力 学 12 111 3 2 12 222 1122 212121 000 000 n n n nn xxx R xxx xxxxxx W WWW WWW WWWWWW （372） 122 2 12 000 000 nn n R W WWW WWW （373） 弱解形式（362）可写成 q TTTT r TdNbdNN UdNqd AAAAA WWWW （374） 由本够关系式（367） ，则进一步有 q TTTT r CUdNbdNN UdNqd A A AAAAA WWWW （375） 由于待定参数与空间坐标无关，可以提到积分号外，因此式（375）可写为U q TTTT r BC BdUNbdNqdNNdU A AAAAAA WWWW （376） 于是，所求的线性方程组为 r K UffM U AA （377） 其中 （刚度矩阵） 22 Tnn KBC BdR A A W （378） （质量矩阵） 22nn MNdR （379） （体力向量） 2Tn fNbdR A W （380） （面力向量） 2 q Tn r r fNqdR A W （381） 方程式（377）是一组线性常微分方程，给定初始条件，可用第 7 章 72 节叙述的数值 积分方法求解。对于静力学问题，因此方程式（377）为一组线性代数方程0u r K Uff A （382 ） 求解线性代数方程组的基本算法参见第 7 章 75 节。不论是静力学问题还是动力学问题， 求解出待定参数后，位移项即可由式（350）确定。u 方程式（377）涉及整个区域，为使微分方程的近似解有满意的精度，人们必须 适当地选取权函数和形函数。对复杂的几何形体，它们的确定通常是一项困难的NWN 工作。在有限元法中，区域有许多简单的微小单元，即有限元来近似。其优点是对简单 的有限单元权函数和形函数很容易确定，然后通过合理的装配，偏微分方程在整个NWN 区域上的解则能够近似地计算。权函数的不同选择导致不同的计算方法，本章NW 33 节简要地讨论了权函数的若干选择。权函数的选择取决于插值的精度、单元的几何形 状以及待定系数的物理意义等因素，这些将在本章 36 节中讨论。待定参数的物理意义取 决于应用的单元类型。对线弹性问题，用有限元法求解位移场时，这些待定参数为结点的 位移。值得注意的是，到目前为止边界条件中只涉及力的边界条件；位移边界0T nqA 条件只有在求解方程式（377）中才被考虑（参见本章节）0uu 本节中方程式由加权余量法导出，更详细的描述可参看Burnett 1987。人们亦可利 用变分法，例如瑞利-迦辽金法来推导有限元方程，如在参考书Reddy 1993中那样。这种 方法由于类似于结构力学中的极值原理，通常受到人们的偏爱，但在应用中变分法对偏微 分方程的要求甚严，且不易被理解和应用。 33 加权余量法 为使方程组式（382）中的待定系数能够唯一确定，权函数必须满足某些要求，u j W 现 代 接 触 力 学 诸如线性独立性。在文献中，人们引进了各种不同的权函数，参见Finlayson Scriven 1966a。本书采用迦辽金法，取权函数为形函数，即 , jj NNNBB WW W （383） 其优点是质量矩阵和刚度矩阵 K 为对称矩阵。在某些情况下，迦辽金法和瑞利-里兹法M 导致同样的结果，参见Reddy 1993。 配点法的基本思想是微分方程只要求在若干点上得到满足。权函数相当于取狄拉克 j W 函数。该方法强迫余量在区域个点上等于零。配点法的(),() jjj xxxx即Wn j x 精度与所选的配点位置密切相关。一般来说配点法中的刚度矩阵和质量矩阵为非 j xKM 对称矩阵。配点法的一种改进方法是子域法，这种方法强迫余量在子域上的积分为n j 零。权函数在子域内为 1，在子域外为零。 j W j j 最小二乘法的实质是使得余量平方和取最小值。它可由变分原理而导出，参见 Papadopoulos 1996a。其权函数。最小二乘法中的刚度矩阵和质量矩阵 j j R u WK 同样为对称矩阵。M 34 基本边界条件 对位移边界条件的处理在Rurnett 1987一书中有着详细的描述，这里仅讨论有关算 法和程序实现的若干细节。位移边界条件可以采用以下方法引入。 1） 矩阵分块法 在用有限元法求解位移场中，位移边界上位移已知，力的边界上力已知，这样 u q 在方程 K UfA （384） 中，每一行要么结点位移已知，要么结点力已知，为了求解未知的变量，方程式 i u i f （384）需要进行适当变动。按结点，位移是已知还是待定需要重新组合方程。 通过调换矩阵和列向量的行以及矩阵的行和列，方程式（384）可用分块矩UfK 阵形式表示，即 qqquqq uquuuu KKUf KKUf A （385） 所有在位移边界上给定的位移变量包含在中，未知的位移变量包含在中，所有在力的边界 上给定的力包含在中，未知的反力包含在中。利用分块矩阵方程的第一块，由下式 qqqquuq KUKUfAA （386 ） qqqqquu KUfKUAA （387 ） 1 () qqqququ UKfKU AA （388 ） 可求得未知的结点位移。 q U 在数值求解中，矩阵的逆无需直接计算，而是采用在第 7 章 75 节描述的方法求 qq K 解线性方程组式（387） 。在解出结点位移后，位移边界上的约束反力可由简单的 q U u f 乘法运算而确定，即 11 () uuqquuuuqqqquqqqquuuu fKUKUKKfKK KKU AAAAAA （389） 在许多有限元程序中刚度矩阵并未建立，实际使用的只是子矩阵，作用力则K qq K q f 通过和作用应的修正。子矩阵和只是在需要计算约束力时才建立。矩阵分 qu K u U uq K uu K 块法虽然原理简单，但其数据结构的操作却相当费时，因为结点位移的顺序性被破坏。 2） 修改元素法 现 代 接 触 力 学 在修改元素法中，结点位移的顺序无需改变，但矩阵和力向量中的元素须作相应Kf 的修改。这里简要地叙述其操作过程：对方程组（384）进行处理。首先检验在第 行中i 是否结点位移已给定，如是，则对力向量中的每个元素分别去掉力的修正项。 i U i f iji K U 在符号上所有未知的结点力则写成。对每一行进行上述操作后，则形成非对称矩阵 i f i U ，修正的力向量和未知变向量。线性方程组KfU K UfA （390） 则可直接求解。方程组的解既包括未知的解点位移又包括未知的结点约束反力。U 现借助一个简单的例子说明这一方法。原线性方程组为 11 22 33 44 15020 90310 242142 40021 K Uf Uf Uf Uf A （391） 1 2 3 4 50201 90310 242142 00214 f U U f A 1 2 3 4 10205.15 031099 0140242.12.414 00211.44 Uf K f U U f A 运用修改元素法，其解则为 1234 11,2,3,10fuuf 修改元素法尽管看上去完美，但由于矩阵的非对称性，在计算中并不比其他方法有K 效，因而很少使用。 在上述两种引入位移边界条件的方法中，方程组需要重新组合，然而可以获得在位移 边界上准确的约束反力。人们亦可采用其他途径，由插值函数的导数和结点位移获得约束 反力的近似解，参见Reddy 1993。这种方法不用重新组合方程组，因此能够节省计算时 间，但以近似误差为代价。如果在时间积分中需要约束反力，则最好运用准确计算约束反 力的方法。 35 坐标变换 通常坐标轴的方向选取与结点位移平行的方向。对图 36（a）中所示的在水平面上 滑动的结点，其边界条件为：待求，已知，待求。P x u x f y f ( )b ( )a P P 图图 3.6 位移边界条件位移边界条件 如果平面如图 36（b）所示，倾斜角为，则自由和约束的位移方向不再平行于坐 标轴。引进新的坐标系，则待求的量为和，已知的量为和。为表达这种坐标 x u y f y u x f 变换，方程式（384）中与点有关的行和列须乘以变换矩阵P cossin sincos （392） 其他所有的行和列则保持不变。借助于正交变换矩阵 1 1 cossin sincos 1 1 P （393） 现 代 接 触 力 学 方程式（384）在形式上可转换为 , TTT KPK P UP U fPfK UfA AAAA （394） 由于矩阵的稀疏性，在上述表达式的数值计算中，自然用不着进行矩阵乘法运算，而只P 需对有关的行和列作相应的修正。 在坐标系中算出的量及通过变换式及即可转换为在,x y U f UP UAfP f A 坐标系中相应的量。通过合理的选择变换矩阵，人们能够方便地处理位移变量间, x yP 的线性相关性，即 iijji aUa Ua UUg kk （395） 一般来说变换矩阵 1 1j 1 1 1 1 1 j iii aa iP aaa 第行 第行 第行 k k （396） 不再具有对称性和正交性。比如，人们要求结点和具有相同的位移，则可转换为 0 iji UU U （397） 在新的坐标系中，则要求位移边界条件必须得到满足，得到解后，通过变换0 i U j U 式则有UP UA 1 jj jjj UPUUA （398） 110 iijji iijijj Up UP UUUUUAA （399） 可见，结点和总有相同的位移。对有些边界条件，比如，人们可以称对地, ij UUU k 处理，例如要求。然而引进这些约束方程时，人们00 ijii UUUUUU和 kk 必须注意约束方程间的独立性和相容性。 36 单元格式 根据式（350） ，位移场由插值函数近似描述。由此产生的误差则取决于单元的选U 择以及插值函数的类型和阶次。下面将对不同类型的单元作一简要描述。 361 等参单元 依据文献Brenner Scott 1994，一个有限单元由三个要素来定义。如果： （1）为一个分段光滑的有界区域（单元区域） ；K R （2）为在域上的一有限维函数空间（插值函数空间） ；PK （3）为函数空间的基（结点变量） 。 12 ,NN NN k P 则称为一有限单元。, ,K P N 在有限元中经常使用等参单元。所谓等参单元是指单元几何形状 xNxA （3100） 和单元的位移场函数 uN UA （3101） 采用相同的插值函数来近似描述，其选择虽然十分自然，但并非必须如此。在非等参单元 中，人们可以对单元几何形状和位移场分别采用不同的插值函数来近似描述。 362 分类 依据文献Dhatt Touzot 1983，有限单元有下述不同的特征： （1）形状（比如三角形单元，四边形单元） ； （2）个插值点的自然坐标；n n （3）自由度个数； （4）结点坐标的定义（拉格朗日单元：位移，Hermite 单元：位移和斜率） ； （5）多项式插值函数的基（比如线性，二次，完整，非完整） ； （6）单元在结点处的连续性（比如， ）符号表示所有至阶连续可微的 01 ,CC s CS 函数空间。 表 31 列出了 5 种平面单元，它们皆为拉格朗日单元，在结点具有连续性。 0 C 表 31 若干单元描述 单元 结点数目 自由度 多项式 基 三角形 3 6线性 1, , 现 代 接 触 力 学 三角形612二次 22 1, , , 四边形48双线性 1, , | 四边形816二次 2222 1, , , | 四边形918二次 22222 1, , , | 363 线性三角形单元 对大多数单元来说，一种行之有效的方法是建立参考单元，然后将真实单元和参考单元通 过一对一对的映射联系起来。这种方法不仅直观而且还可以避免不必要的重复运算。图 37 演示在参考区域平面内三角形单元和它在真实区域平面内的映射。参考单元由 r V 不等式确定。通过映射01,01,1 123 ():( )( ;,) re VVxxxx x x （3102） 可确定参考单元相对应的真实单元。 e V 人们亦可把坐标系看作为描述真是单元的非直角局部自然坐标系。对三角形, e V 单元，由式（3100）和（3101）可得 1 1232 3 ( ), x xNxNI NI NIxuNU x AAA （3103） 其中线性插值函数需满足下述条件： N 1 0 ij ij N ij 当 当 （3104） 表 32 列出了这些插值函数及它们对自然坐标的导数。 表 32 线性三角形单元插值函数及其导数 / 1111 210 301 iii iNNN 三角形单元的插值函数被描绘在图 38 中，他们具有下述重要特征： （1）参考单元的结点被映射到真实的结点 1 21 3 10000 x xIxx x A结点：， 1 22 3 1000 x xIxx x A结点2：， 1 23 3 0100 x xIxx x A结点3：，