第一章算法 §1二维 Laplace方程外问题 设r是平面上的一个凸多边形,它的外部区域是Ω(图1)。 我们在Ω上求解 Laplace方程的 Dirichlet问题 △=Q (x,y)∈D,(1.1 (1.2) 或者将边界条件(1.2)换成 01 (1.3) 01 (I.1),(1.3)是…个Nem 题,以后我们永远以ν表示外法 线方向。与(1.1)联系的应变能 是 dxdy (I.4 为求解问题(1.1),(1.2)或者(1.I),(1.3),我们将2剖分为 无限多个三角形单元,不妨设巫标原点O在F内部,取一个常 数5>1,以(点为相似中心,以§,52,…,5k,…为比例常数,作 F的相似形,它们分别记作,F2,…,Ik,…,每两个多边形之 间的ⅸ城称为一“层”(图2)。然后,我们将每一层进一步剖分 成单元,可以采用如下方式:在!上取若干点作为节点,其 中顶点必颁坛节点,在毎条边上还可以按需要适当地再选一些节 点,点O点出发作射线与各节点相连,这样就将每→层削分成相 似的若个四边形,再将每个四边形潮分成两个三角形,需要注意
约是,每层的剖分方式必须一致(图2)。 边界r的形状可能不是一个凸多边形,而是非常复杂的形 状。这时我们可以将厂。的外部区域分解成一个形状规则的无界 区城与一个形状复杂的有界区域。无界区域仍记作9,它可以按 照上述办法称剖分。对于有界区 域,则可以按照常规的办法作有 限单元剖分(图3)。我们还使这 两种剖分互相协调一致,即在交 界处,三角形单元的节点与边互 相吻合。在这种情况下,我们仍 然先分析区域Ω,我们将给出2 上的组合刚度矩阵Kz的箅法 Kz是一个有限阶的矩阵,在使用的时候,或者用它直接求解边 值问题,或者将它与其它单元的刚度矩阵叠加,以求解形状复杂 的区蚁的边值题。对于后者,Ω是当作一个单元来处理的。 现在考虑多边形fk。从某一点开始,按逆时针方向,将Fk
上的各节点排一个次序。设节点共有n个。(1.1)的解在节点 上的值毡排了一个次序,记作y4),k2),…,3),我们将它们排 列成…个列向量(v41,y2,…,y)“,记作Jk,这里T表示矩阵 或者向量的转置。 考虑多边形Fk1与k之间的第七层。对于每个三角形单 元,在节点上的值给定以后,作线性插值,可以按常规的方法求 单元刚度矩阼(例如参看52]),将第k层上的各单元刚度矩阵按 节点叠加,可以得到一层上的斓度矩仵。我们将这个刚度矩阵记 K (15) 即如果第k层上的应变能是Wk,则 K A (v-1, K 其山K…,K′,A都是n阶矩阵。从有限元方法的基不理论可知 (1.5)遠对称知阵、国此1与A互为转,K,K拊是对称 从(1、)和剐分的似性,不难判断、各层的飓度姸阵是札 的。将犴层的刚矩阵按节点叠加,可以得到一个无穷阶的总 树度絆阵。考虑到方程与边界条件,可以给出一个无穷阶的代数 程组。令K=K+K0、则有 Kayo -Ay,=Jos (⊥.6 A,ye-Ky1-A (1.6冫: 要鲁看BP咖D非t。。看面 Al K y (1.6) 物e专·D·卷·些脂 当考虑问题(1.1),(I.2)时,在F上不要列出方程,因此代要 (]6)。·而y是已知的。当考虑问题(1.1),(13)时,在上
要列出方程,即方程(1.6),其中∫是“等效节点力”,用常 规的有限元方法,它是可以从边界值∫求出来的。目前,我们的 目的是给出口上组合刚度矩阵的算法,因此可以暂且不管上述 边界条件的差别。 我们将在下章证明,存在一个n阶实矩阵X,使 Uk+1=Xy (1.7) X称为转移矩阵。以(1.)式代人(1.6)得 K-A+KX-AX)yo=0 但是y0可以是任意向量,所以X满足方程 A TX2-KX+A-0 (1.8) 令 K A R A 0 RI 其中是单位阵,方程(1、6)k可以写成如下形式: J R R (1,10) 设A,9分别是矩阵X的特征值和特征向量,则有 X9=Ag (1.1L) 在(1.10)中取无=1,令孙=9,并且注意到(].7).(.I1)得 g R1 nR (1.12) g 因此)与j 分别恳矩阵束R1-λR2的广义特征值与广义特 征向量1},为了求矩阵￥,我们先解上述广义特征值问题,求 出它的全部2n个特征值与特征向量,然后丛中选出ⅹ的n个特 征值与特征向量
以后,我们总是以det装示一个矩阵的行列式。我们计算 det(r-2R2)=det /x-2.AT A λ 在以上分块矩阵中,以λ乘第一列并加到第二列上得 K -A+AK-ALAT dat(R1-λR2)=d (-1)det(-A+1K-22AT) 由K的对称性,上式是一个λ的对称多项式,因此,如果九≠0是 个特征侑,1/也是一个特征值。我们将在下章中证明:X的 特祇值满足I,而且当14}=1时,必然有=【,其对应的 特征向量9;=(1,1,…,1)T,因此,我们只需将矩阵束R,-λR2的 特征值中满足!∴<I的选出,再添上λ=1、设特征佤为λ1,λ ·,九n,对巨的特征向量是只:,92…,9n,令 T=(;,9 A=diag(λ:,λ2,…,λn), 其中d表示对角矩阵或者对角矩阵。由(1.11)得 XT=TA Y: TAT (1.13 (1.13)就是计算转移矩阼X的公式。以(1.7代入(1.6)得 Axo=f 1.14) 其中Kz=A-AX,它称为组合刚度矩阵,我们将在下一章中 证明Kz是一个对称半正定矩阵 红合刚度姸阵与总刚度矩阼是有区别的。在现在的情况下, 总刚度矩阼是个无穷阶矩阼,通过它可以得到任意节点值产生 的应变能,而組合刚度矩阵仅是一个n阶矩阵,通过它只能得到 当节点值y,y1,…满足方程(I,6)4,(1.6)2,…时的应变能。从总 刚度矩阵到组合谢度矩阵的过程可以看作是一个消元过程。通过 这个过阻 都消去了 利用无限元方法求解的过程是这样的;先求出转移矩阵与组 5