（计算数学专业论文）用边界元法求解三维摩擦型弹性接触问题.pdf

重庆大学硕士学位论文 中文摘要 中文摘要 本文提出了一种基于迭代技术的边界元法求解三维无滑动摩擦型弹性接触问 题。利用广义格林公式和基本解得出弹性问题的边界积分方程，并通过选取适当的 局部坐标系来表达接触条件以便于将不同弹性接触体的边界积分方程祸合。采用迭 代法，通过寻求与接触条件相协调的接触边界位移及面力增量来确定无滑动接触区 域的大小。由于摩擦的存在，考虑到加载历程的不可逆性，本文采用模拟外部荷载增量 的方法来追踪加载历程。数值计算采用配置法计算耦合的边界积分方程，使用四边 形常单元。奇异积分用解析积分公式计算，文中推导了相应的积分公式。 论文采用f o r t m h 9 0 语言编制了三维无滑动摩擦型弹性接触问题的通用边界元 程序，并对相容和扩展型接触问题进行了数值计算。数值实验的结果表明文中所提 算法是收敛且是行之有效的。 关键词：边界元法接触问题迭代法增量理论 重庆大学硕士学位论文 英文摘要 a b s t r a c t ab o u n d a r ye l e m e n tm e t h o db a s e do na l li t e r a t i v et e c l m i q u ef o r s o l v i n gt h e t h r e e - d i m e n s i o n a le l a s t i cc o n t a c tp r o b l e mw i t h o u ts l i d i n gi s p r e s e n t e d t h eb o u n d a r y i n t e g r a le q u a t i o nf o re l a s t i c i t yi sd e r i v e dt h r o u g ht h eg e n e r a lg r e e n sf o m m l aa n dt h e c o r r e s p o n d i n gf u n d a m e n t a ls o l u t i o n t h ep a p e rr e p r e s e n t st h ec o n t a c tc o n d i t i o n s ，w h i c h a r ee s s e n t i a lf o rt h ec o u p l i n go ft h eb o u n & l i - yi n t e g r a le q u a t i o n so ft h et w od i f f e r e n t e l a s t i cc o n t a c tb o d i e s ，i nal o c a lc o o r d i n a t es y s t e mp r o p e r l yc h o s e n i nt h i sp a p e r , a l l i t e r a t i v ep r o c e d u r e ，w h i c hs e e k st h ei n c r e m e n t so f b o u n d a r yv a l u e sf o rd i s # a c e m e n t sa n d t r a c t i o n sc o m p a t i b l ew i t ht h ec o n t a c tc o n d i t i o n s i sa p p l i e dt od e t e r m i n et h ec o r r e c t c o n t a c tz o n ew i t h o u ts l i d i n g d u et ot h ee x i s t e n c eo ff r i c t i o n ，t h ep a p e ru t i l i z e sam e t h o d o fs i m u l a t i n gt h ei n c r e m e n t so ft h ee x t e r n a ll o a dt o 臼粼t h eh i s t o r yo fl o a d i n gt a k i n g a c c o u n to ft h ei r r e v e r s i b l ec h a r a c t e ro ft h eh i s t o r yo fl o a d i n g i nn u m e r i c a lc a l c u l a t i o n , t h ep a p e re m p l o y st h ec o l l o c a t i o nm e t h o dt oc o m p u t et h ec o u p l e db o u n d a r yi n t e g r a l e q u a t i o n sw i t hq u a d r a n g l ec o n s t a n td e m e n t s t h es i n g u l a ri n t e g r a l sa r ec a r r i e do u tb y a n a l y t i c a li n t e g r a t i o n ag e n e r a lf o r t r a n 9 0p r o g r a mf o rc o m p u t i n gt h e3 - df r i c t i o n a le l a s t i cc o n t a c t p r o b l e mw i t h o u ts l i d i n gb yb o u n d a r ye l e m e n tm e t h o di sc o m p i l e da n ds o m en u m e r i c a l e x a m p l e si n c l u d i n gc o n f o r m i n ga n da d v a n c i n gc o n t a c tp r o b l e m sa l eo f f e r e dt oe x a m i n e t h ec o n v e r g e n c ea n df e a s i b i l i t yo ft h ed e s c r i b e dm e t h o d t h en u m e r i c a le x p e r i m e n t s i l l u s t r a t et h a tt h em e t h o dp r e s e n t e di nt h ep a p e rw o r k sw e l l k e y w o r d s ：b o u n d a r y e l e m e n tm e t h o d ；c o n t a c tp r o b l e m ；i t e r a t i v em e t h o d ； i n c r e m e n t a lt h e o r y i i 重庆大学硕士学位论文1 绪论 1 绪论 数学物理问题，如油、气藏的勘探与开发，大型结构工程，航天器的设计，天 气预报，反应堆等涉及到传热、传质、渗流等问题的计算，无不归结为求解大型偏 微分方程，当区域大且形态极不规则时，求问题的解析解是根本不可能的，所以数 值计算成为求解此类问题的有效方法。边界元法是随着有限元计算方法的发展而在 七十年代后期被提出来的求解偏微分方程初边值问题的数值计算方法，由于它一般 只需进行边界剖分，有降低一维进行数值计算的特性，所占计算机内存小，计算时 间省，与有限元法相耦合能很好地解决工程实际问题等优点。目前，这一方法已在 各个工程领域如：力学、热学、航空动力学、机械设计与加工等获得广泛应用，特 别是边界元法中的基本解性质往往能恰当地反应出无穷远处的边界条件，使得利用 它来研究无限域问题具有其它数值方法所不能具有的一些独特的优点，尤其近十年 来边界元法在研究非线性和各向异性问题方面取得了很大进展。使得边界元法在工 程中很受青睐。 所谓弹性接触问题，即求解两个或多个与接触条件耦联的n a v i e r 方程，是弹性 理论中一种常见的且非常重要的一类问题。自从1 8 8 1 年h e r t z 发表了关于弹性体之 问法向接触问题的著名论文之后，在该领域，科学家展开了大量的研究，包括理论 和实验方面。由于在实际工程中，接触体的几何形状以及所受荷载等均很复杂，这 时寻求解析解几乎是不可能的。因此必须采用数值解法对该问题进行数值求解。当 前使用较为广泛的数值解法主要有有限元法，边界元法等。 1 1 研究现状评述 1 1 1 边界元法发展现状 工程，物理问题的数学模型的解如果可以转化为某个边界积分方程的形式，那 么用数值方法求解边界积分方程以达到求解数学模型的方法就是边界元法，所以边 界元法有时被称为边界积分方程法。它是将区域内的偏微分方程的边值归化为边界 上的积分然后在边界上离散化求解的一种数值方法。其基础在于边界归化，如果 边界归化的途径不同，可以从同一边值问题得到几个不同形式的边界积分方程，所 以可以导致不同的边界元法，目前大约有这几种：采用g r e e n 公式归化的直接边界 元法( 本文即采用此方法) ，采用位势理论归化的间接边界元法，还有用g r e e n 函数 代替基本解的自然边界归化的则称为自然边界元法。 间接归化法的提出，可以追溯到1 9 6 3 年j a s w o n 的论文，他首先对边界元的间 接法提出了完整的概念，间接归化是从基本解及位势理论出发得到f r e d h o l m 积分方 重庆大学硕士学位论文 l 绪论 程，这种方法所求的物理量不是原问题的解，而是根据线性叠加原理在边界上分布 的虚拟的“源”分布密度，它本身只有间接的物理意义，在位势理论中可用单层位 势或双层位势，甚至用这两种位势的组合来表示所求的物理量。 直接边界归化法则是从基本解和广义g r e e n 公式出发，将微分方程的边值问题 转化为边界上的积分方程，工程界常用的所谓，“加权余量法”就属于这一类型，与 间接法不同的是，直接法不引入新的变量，积分方程的解就是原问题未知量的边界 值或与边界上未知量的法向导数相关的值，所以这一方法在使用上比较方便丑易于 理解，故在工程界很受欢迎。 我国学者冯康1 9 7 5 年提出了利用g r e e n 函数从g r e e n 公式出发将偏微分方程的 边值问题转化为边界上的含有发散积分的有限部分的积分方程的思想，这就是正则 边界归化，或称自然边界归化。处理利用这种转化方式而得到的含强奇异积分方程 的数值方法称为正则边界元法，后来称为自然边界元法。这种方法与有限元方法的 耦合，开辟了处理实际问题的前景并弥补了单独用有限元法难于处理无限区域问题 的缺陷。 早期的边界积分法只注重把偏微分方程的边值问题怎样归化为边界积分方程。 近代的边界元法则在把偏微分方程的边值问题归化为边界积分方程的基础上引进了 有限元的离散技术，着重于积分方程的数值求解过程。英国s o u t h a m p t o n 大学的 b r e b b i a 教授1 9 7 8 年提出边界元方法，( 事实上b e m 是美国的t a c r u s e 教授于 1 9 7 2 年首次提出，只是后来他又放弃了这种称法改称其为边界积分方程法) 表明边 界元方法作为一种与有限差分法，有限单元法并列的求解偏微分方程的数值解法已 取得了独立的地位。 边界元法的提出离不开计算机的发展，在计算机未得到普及之前，对边界积分 方程的研究只限于理论研究。二十世纪四五十年代，随着电子计算机的出现，发展， 大规模数值计算成为可能。使得数值计算得到了前所未有的发展。可以说数值计算 的发展与计算机的发展是并驾齐驱的。甚至可以说，没有计算机就不可能有数值计 算的今天，也不可能产生新的计算方法。 二十世纪四五十年代出现的有限元法作为一种先于边界元法的一种数值方法， 为边界元的出现和发展奠定了基础。有限元法是基于变分原理，采用分片插值函数 逼近，先把整个区域进行单元剖分后再进行计算。正是因为计算机和有限元法的出 现。大量的科学研究可以通过数值分析和数值模拟来进行。到了七八十年代，其理 论已趋于成熟，其应用也遍及科学，工程的各个领域。 边界元在有限元比较成熟的时期出现，得到迅猛发展，并在某些方面可取代有 限元。与有限元法相比它有如下优点： 1 ) 有限元必须在整个区域上进行单元剖分。边界元一般只在边界上进行单元划 重庆大学硕士学位论文 1 绪论 分，这可降低一维进行计算，因而在数据准备上边界元法要少很多，计算量也大为 减少 2 ) 工程中的无限域或半无限域问题，用有限元来处理，则是截取一个较大的计 算范围，在外边界给出近似的边界条件。为了保证计算精度，范围越大越好，但是 计算范围取大，结果是单元增多，存储容量和计算时间都相应增加，八十年代初期 出现了一种无限元或无界元的方法，虽可以部分解决这类问题。但仍要把计算范围 取得足够大，然后在外围再布置一层无限元，所以工作量减少是有限的，而边界元 法中由于采用的基本解本身就适用于无限域或半无限域，因此求解此类问题时不必 再需要“外边界”，从而使问题大为简化。 3 ) 工程中常有所谓的奇异性问题，如裂纹尖端的应力，应力集中，有源的渗流 场问题，用有限元法来求解这类问题时，必须将奇点附近的网格划得很细。因而数 据信息和计算量都很大。例如坝基和坝体有若干排水孔，而其尺寸比结构尺寸小得 多，用有限元法来计算时，为了真实反映排水孔的作用，必须将孔附近的区域划分 为很细的网格，但又不能将整个区域分得很细，否则存储量太大以致不能计算，由 孔附近较密的网格变到边界附近较稀的网格，需要很多过渡单元，但是自动划分技 术却无能为力，所以数据准备工作量就大大增加，如果用边界元来处理这一问题， 由于排水孔的边界和结构的边界剖分彼此互不干扰，不要求用过渡单元来协调，问 题处理起来就很简单了。 4 ) 流体力学中的自由面问题，由于自由面的位置常难以预先估计，只能任意假 定初始面经过选代计算来求解，用有限元法来求解不仅程序复杂，而且在迭代过程 中如果自由面过大变动会使一部分单元随着有过大变形，影响计算。边界元方法因 仅要求对边界进行剖分，迭代过程中自由面变动的程序处理就很简单。 5 ) 由于边界积分方程本身是所求问题的一种精确表示法，其误差仅来自边界 离散化的处理，数值上的近似计算也只是因为在边界上无法闭合积分解而引起的， 若采用相当成熟的数值积分法，比如采用曲边单元和在边界上连续变化的函数分布， 则积分引起的误差可以很小。采用边界元法，一旦求得边界值，则可以由积分表达 式解析地求出域内解，且处处连续，因此精度很高，这是比有限元法好的地方。 正是由于边界元法具有以上的优点，所以备受工程师们的喜爱，使得边界元法 得到迅猛发展，自1 9 7 8 年第一届国际边界元法会议召开以来，每年都要召开一次， 其研究领域已遍及抗体力学，弹性力学，弹塑性力学，岩土力学，热学，声学，电 磁学，生物细胞学，航空动力学等。研究的问题也越来越复杂。早期的边界元法研 究主要是线性问题。而目前主要趋于以下几个方面。 a ) 强非线性问题 边界元法的一个特征是利用基本解得到积分方程然后求解。对线性问题一般都 重庆大学硕士学位论文 i 绪论 能得到基本解，对于非线性问题难于得到用于数值计算的基本解，但仍有办法处理， w u ，b r e b b i a , b r u c h ，l i g g e 等曾做过很多工作，取得了一些进展。但仍有很大的研究余 地。 b ) 变系数问题 系数是与坐标有关或与坐标和时间都有关的问题，一直是边界元法应用的弱点， 因为此类问题同样难以求得可用于数值计算的基本解。目前的处理方法是采用将系 数矩阵通过特征求值，坐标变换或近似基本解等方法将问题转化为常系数或方向性 问题，使研究大为简化。 c ) 多区域问题。 在工程中有些问题的各个区域的物理性质不相同，从而导致区域分解算法或区 域耦合法。将一个区域分为若干个子区域联立来求解，虽然增加了边界单元的数目， 但各个区域仅仅在交界上相互联系，其线性代数方程组的系数矩阵是带状稀疏的， 也就是可以减少计算机存储量。这种方法一般适应于多介质问题，形状不规则 的区域。 对于多区域问题用边界元法求解，已有不少研究，除了完全采用边界元技术外， 还有边界元法和有限元法或其它计算方法的祸合方法。如要在边界元法中引进并行 算法，则要采用区域分解算法，区域分解法不同于区域耦合法，它是在每个区域上 单独求解，正是边界元与有限元耦合法或其它方法和边界元结合法拓展了边界元法 的应用范围。 总之，边界元方法自七十年代末提出以来，在二十几年内取得了很大进展，在 计算方法领域占有很重要的地位，许多关于边界元法的专著己出现。边界元法作为 一种独立的计算方法在工程计算中得到了广泛应用。 1 1 2 弹性接触问题的数值解法发展现状 在土木，建筑，水利工程，石油化工，机械工程等领域中广泛存在着接触问题。 例如。基础梁与地基之间：混凝土大坝出于施工和散热的需要而设置的临时或永久 性的分缝：齿轮传动中的轮齿接触等。对于一些几何形状规则的物体间的接触问题 ( 如两个弹性球体之间的接触问题) ，能用经典的数学工具( 如积分方程法) 来求解。 虽然经典接触力学得到的解的结果很精确，但是其解决问题的范围却十分有限， 在实际工程中，对接触问题的求解大多采用数值解法，有限元法与边界元法是两种 当前被最广泛地用来求解弹性接触问题的数值方法。 1 ) 有限元法在接触问题中的应用始于6 0 年代术。1 9 7 0 年，w i l s o n 和p a r s o n 首先研究了二维弹性无摩擦接触问题的有限元解法。继而c h a n ，t u b a 以及o h t e 等人先后将有限元法推广到带摩擦的二维和轴对称的弹性接触问题。但上述工作均 未考虑加载过程中的不可逆性。在工程接触问题中，通常存在有摩擦力，而摩擦力 4 重庆大学硕士学位论文 1 绪论 所做的功一般体现为能量耗散，因此若不考虑摩擦力的影响或不考虑加载过程中摩 擦力的不可逆性，就难以准确的反映工程问题( 同样也是边界元法求解接触问题所 需要解决的) 。1 9 7 3 年，t s u t a 等人提出了一种基于载荷增量理论的有限元法，用于 求解带摩擦的接触问题，较好地解决了加载过程中的不可逆性，并成功地作了二维 算例分析。1 9 7 6 年f r e d r i k s s o n 等人从理论上进行了较严格的推导，建立了弹性接触 体的增量控制方程，并用有限元位移法进行求解。1 9 7 9 年o k a m o t o 和n a k a z a w a 等 人从力学基本变分原理( 虚功原理) 出发，建立了增量控制方程及有限元解法。 七十年代末期以后，基于有限元法的一些数值方法相继出现，具有代表性的主 要有以下几种方法： a ) 传递矩阵法。该方法是一种直接法，它直接根据接触条件通过迭代来决定 节点的接触状态； b ) 间隙有限单元法；它在接触体系的间隙之间引入一种人为添加的，虚拟的， 只能受压；不能受拉的间隙单元。它的开和关取决于接触物体的特性和载荷情况， 通过迭代调整该单元的压缩弹性模量。间隙单元组成的子区域对应的矩阵方程经过 变换后和其它单元的矩阵方程具有相同的形式。此方法概念简单，便于使用，但在 处理复杂的实际问题时收敛性不能得到保证； c ) 罚函数法：也称为罚有限元法。该方法取符合一定条件的接触单元为罚单 元，将接触区的非嵌入条件作为罚项引入接触系统的总势能中，将约束变分问题转 化为罚优化问题来求解。罚单元的节点可以由属于不同单元的节点组成，也可以由 某个实际单元内的一些节点之间的约束来构成。因为这个原因，罚单元在接触问题 中应用起来有其优越性。但是对罚因子的选取直接决定了方程的求解精度，因此有 对会造成病态方程而无法求解。 d ) l a g r a n g e 乘子法；此方法将接触问题考虑为有约束的最小值问题，通过引 进l a g r a n g e 乘子来构造l a g r a n g e 函数，使上述问题转化为l a g r a n g e 函数的无条件 极值问题；将一组变分不等式化为适合有限元计算的一组等式来对接触问题求解。 除了上述基于有限元法用于求解接触问题的数值方法之外，还有一些有效的数值方 法，如：基于柔度法的混合法( 也称为广义子结构法) ；在l a g r a n g e 乘子法基础上 发展起来的摄动l a g r a n g e 法和增广l a g r a n g e 法；以及数学规划一有限元法等。 总之，当前利用有限元法求解弹性接触问题不仅在数值方法上而且在数学理论 上已经比较成熟，并且，由于有a n s i s 等基于有限元分析的可对实际工程问题进行 数值模拟的大型软件包的存在，有限元法在工程中得到了广泛的应用。 2 ) 对于接触问题，边界元法相对有限元法有着明显的优点。使用边界元法，问 题被有效地降低一维，求解问题只需对边界进行离散化处理，对于绝大多数应力问 题最大应力出现在边界上这一点显得尤其有效，因为此时，不需要计算内点的应力 重庆大学硕士学伉论文 1 绪论 和位移值，从而在很大程度上减少了计算时间，降低了问题的复杂性。有限元法不 能直接利用接触条件，这就意味着接触压力不能作为一个独立可求解的变量出现在 求解公式中，而在边界元法中与接触区有关的变量，如面力、位移，不仅可以在边 界积分方程中得到直接处理，而且在数值上可以达到相同的精度。另外边界元法处 理应力集中问题更为有效，并可求解任意一点的应力等参数。虽然在处理非线性问 题上，边界元法有某些缺陷，如偏微分方程不能完全精确地被转化为边界上的积分 方程，材料的非线性性等。但是特殊地对于线弹性接触问题，边界元法仍然是有效 的，因为此时非线性仅出现在边界上。近年来，对于接触问题边界元法的研究越来 越受到学术界与工程界的重视。 将边界元法首次应用于接触问题的研究是瑞典人a n d e r s s o n t 【7 l 【b 】。1 9 8 0 年， 他从全量理论出发，提出了二维弹性无摩擦接触问题的边界元法，推导了相应的边 界积分方程，并作了算例分析。但是该工作没有考虑加载过程的不可逆性。次年， 他吸收了有限元技巧，从增量理论出发，导出了边界元增量求解的迭代方程，将边 界元法推广到二维弹性带摩擦接触问题，并考虑了加载过程的不可逆性。继而，针 对二维问题，讨论了接触区面力不连续点的处理问题。1 9 8 2 年，a n d e r s s o n t 采用 线性和抛物线等高次边界元对弹性基础与卷轴及平面冲压机的接触问题进行了数值 分析，与f r e d r i k s s o n l 9 1 采用有限元方法所得结果取得一致。为了避免由使用二次抛 物线单元所产生的几何不协调性，他提出了对边界单元的中点添加线性约束的方法。 1 9 8 5 年，p a r i s f t 5 6 1 采用不连续单元概念处理接触区面力不连续问题( 如应力集 中问题) ，并用增量迭代技术求解。1 9 8 6 年，a d u l m i h s e i n 【5 7 1 使用二次等参元和增 量理论对轴对称问题进行了研究，1 9 8 7 年，k a r a m i t 5 ”采用二次等参元分析了具有 体力的摩擦型热弹性接触问题；j i n t 5 9 】利用o d e n 和p i r e s 提出的非局部以及非线 性摩擦定律对二维摩擦型接触问题开展了研究。8 0 年代末期，p a r i s f i ”j 等人对多体 接触且具有多个外部荷载的接触问题开展了研究，解决了施加全新荷载时所产生的 不协调性。c j a k w a y a a 利用边界元法用n a v i e r - s t o k e s 方程代替r e y n o l d s 方程对热流 体弹性动力接触问题进行了数值求解。a t r a l l i 和c a l e s s a n d r i 从基于接触压力的变 分不等式出发，利用直接边界元法对无摩擦的接触问题进行了数值分析。k u i c h 提 出了将柔度法和边界元法相结合的思想来求解接触问题，并强调了采用不连续单元 的优点：继而，s t a k a h a s h i 和b r e b b i a 采用边界元法来计算柔度矩阵并取得了很好 的数值结果。 以上工作仅涉及到二维弹性接触问题。9 0 年代，对于三维弹性接触问题的边界元 法研究也取得了一些进展。h a s e g a w a 针对刚性轴横向插入空心圆柱体的孔中而产 生扭转的问题，利用g r e e n 公式将扭转体力的区域积分转化为对边界的积分，分析 了轴与孔之间无摩擦的三维弹性接触问题，并进行了轴对称体无摩擦接触问题的边 重庆大学硕十学位论文 1 绪论 界元法研究。l i us h u b i n 等人用全量法推导了无摩擦三维弹性接触问题的边界积分 方程及离散方程，并利用凝聚技术，得到仅包括接触边界上的未知量的迭代方程， 大大有利于提高迭代计算效率及减少计算机内存量。1 9 9 2 年和1 9 9 4 年，西班牙 p a r i s f 帅j _ 【“1 等人利用迭代法和增量法分别对三维无摩擦和有摩擦接触问题作了数 值分析，完整的使用了库仑摩擦定律，利用n e w t o n r a p h s o n 方法解决由于滑动所引 起的非线性性。 3 ) 最后，由于有限元法在处理非线性材料构成的弹性体接触问题以及含有薄板 的接触问题上较边界元法更加有效，科学家们在耦合有限元法和边界元法求解弹性 接触问题方面也开展了研究。a l a n d e n b e r g e r 2 1 1 利用间隙有限单元对弹性接触问题 进行了边界元数值分析；y e z a w a 和n o k a m a t oi 2 0 l 利用边界元法和有限元法研究 了磁头与磁盘之间的接触问题。 1 2 本文使用的方法 在大多数接触问题中，接触区域是外部荷载作用的函数，即接触区域随外部荷 载的变化而变化。当需要考虑摩擦时，位移和摩擦区域的最终状态用加载历程来确 定，因此必须追踪整个加载历程。本文研究的问题是三维无滑动摩擦型弹性接触问 题。利用位移和面力接触条件耦合边界积分方程来解决单个边界积分方程本身所具 有的非线性性质( 因为接触区域的边界函数全为未知函数) 。通过一迭代过程来寻求 与迭代接触条件，即不属于接触区域但距其很近的点对之间不发生相互渗透，属于 接触区域的点对法向应力为压力，相协调的边界解来确定接触区域的大小。对于扩 展型接触问题通常选取比实际接触区域大的区域用于首次迭代，通过迭代运算，用 于下一次迭代的接触区域逐步减小以至最终达到实际接触区域。外部荷载增量的大 小用模拟的方法来确定即按比例指定荷载增量来追踪加载历程，并通过迭代法计算 位移和面力未知量的增量。数值上使用配置法采用四边形常单元进行计算。 本文的主要研究工作有： ( 1 1 详细推导三维弹性接触问题的边界积分方程； 佗) 推导了三维弹性问题奇异积分的解析解法： ( 3 ) 给出了确定接触区域大小的迭代过程： ( 4 ) 给出了弹性接触问题基于迭代法的边界元方法： r 5 1 用f o r t r a n 9 0 语言编制出利用边界元法求解三维无滑动弹性接触问题的通 用计算程序； ( 6 ) 给出了一些具体的算例，证明了算法的实用性和有效性。 重庆大学硕士学位论文 2 三维弹性接触问题边界元解法 2 三维弹性接触问题边界元解法 2 1 问题的提出 图2 1 1 问题的描述 f i g2 1 1t h ed i f i n i t i o no f t h ep r o b l e m 考虑两个三维弹性接触体d ( 丘= 彳，b ) 其边界分别记为s ( k = a ，b ) ( 参见图2 1 1 ) ，坐标系采用笛卡儿坐标系， 每个弹性体的边界按接触情况分为 磷一接触区域 彰一非接触区域( k 利，召) 显然有s = s 5 u 醛( k = a ，b )( 2 1 1 ) 对每个子区域掣，黔又可以分为 s 5 = s 己u s 盖 ( k 刊，口)( 2 1 2 ) s l 一接触粘贴区域( 其中的点对没有发生相对位移) s l 一接触滑动区域( 其中的点对由于外部荷载作用发生了相对位移) 磷；5 品u 芦各( k = 彳，b ) ( 2 1 3 ) s ：，一位移为已知的区域 s s 一面力为己知的区域 假设接触体满足小变形条件，在本文中假设无滑动摩擦， 8 重庆大学硕士学位论文 2 三维弹性接触问题边界元解法 故有 s ：z s g * s c l s c = s l ( k 刊，丑) ： ( 2 1 4 ) ( 2 1 5 ) 昂实际接触区域； 为方便起见，不考虑体力作用。 考虑弹性问题 k ，j = 0 i nd ( 1 ) ( p 1 ) 盯：= c “i 。 i nd ( 2 ) ( k = a ，b ) “? = “?o n 踮 ( 3 ) i j = l ，2 ，3 t = 寸 o n 蹬 ( 4 ) 仃：n j = | ；，破= h ：，i s i o n s ；( 5 ) 其中k 一弹性体d 上点的应力张量，且有盯乞= a 盯； “! 一d 上点的位移向量： f ? 一d 。上点的面力向量5 n 一边界点上的单位外法向量； c ；t m 2 凳6 u 6 h + l k 嫡d 6 i m 七6 t 小 ，t l a m e ( 拉梅) 弹性常数； 最- - k r o n e c k e r 符号； ( 1 ) 一( 5 ) 式引用了逗号指标符号及e i n s t e i n 求和约定。 2 2 一般弹性问题的边界元解法 在本节，首先研究单个弹性体边界元解法并将弹性问题采用算子方程的形式表 不。 考虑问题 ( p 2 ) l u = 0 “i l = “ t ir = i 加d o n f i f = a d f = f 1 + f 2 o n f 2 其中：l = 肚+ + ) g r a d d i v 对问题( p 2 ) 应用广义格林公式和基本解，即可得s i m i g l i a n a 位移公式 “( p ) = i u 口( b q ) f ，( q ) ( 恣( q ) 一亿( 弘q ) u j ( 9 ) 勰( q ) r r p dq f 其中u p ( p q ) ，乃( p ，q ) 为k e l v i n 基本解，其表达式分别为 9 重庆大学硕士学位论文2 二维弹性接触问题边界元解法 故有：s jz 鲜z s ，；( 2 1 4 ) 品= s 占( k = ； ( 2 1 - 5 ) s ，一实际接触区域； 为方便起见，不考虑体力作用。 考虑弹性问题 口! = 0 i nd ( 1 ) ( p o盯：= c “i 。 i nd ( 2 ) ( 世。，b ) “! = “； o n s 品 ( 3 ) i j 2 l ，2 - 3 t = t t o n s 鼻 ( 4 ) 砖”。= 寸，“? = 群，i ，o n 研( 5 ) 其中弹性体d r 上点的应力张量，且有a 嘉= 。； “一d 。上点的位移向量： ，一d 。上点的面力向量； n 一边界点上的单位外法向量； c ：- 。2 凳6 v 6 t 。+ p 够n 6 m + 6 。6 i i 、 “一l a m e ( 拉梅) 弹性常数； d _ k m n e c k e r 符号； ( 1 ) 一( 5 ) 式引用了逗号指标符号及e i n s t e i n 求和约定。 2 2 一般弹性问题的边界元解法 在本节，首先研究单个弹性体边界元解法并将弹性问题采用算于方程的形式表 示。 考虑问题 ( p 2 ) 加d o n lf = o d , f = f 1 + r 2 o hl 其中：三= 肚+ + l t ) g r a d d i v 对问题( p 2 ) 应用广义格林公式和基本解，即可得s i m i g l i a n a 位移公式 “( p ) = p 。( p ，q ) f ，( q ) d s ( q ) 一j 毛( p ，q m ( q ) 勰( q ) rr p d q r ( ，2 - 1 ) 其中u ，( bq ) ，瓦( p ，q ) 为k e l v i n 基本解，其表达式分别为 其中u ，( bq ) ，瓦( p ，9 为k e l v i n 基本解，其表达式分别为 9 嚣一 工”r 重庆大学硕士学位论文2 三维弹性接触问题边界元解法 u 。( b q ) = 而丽1 丽r 3 4 u ) 磊+ 。r ，1 ( 2 2 2 ) t o ( p , q ) 一南毫 ( 1 _ 2 v ) 毛+ 3 分啊( 1 - 2 v ) v ，1 j ( 2 2 3 ) 其中，= i p o l y 泊松比 e 一弹性剪切模量 铲， 令p d 斗pe f ，进一步可得相应于问题( p 2 ) 的边界积分方程 c o ( p ) u i ( p ) + j 毛( p ，q ) u j ( q ) d s ( q ) = j u , j ( p ，q ) f ，( q ) 您( q ) ( 2 2 4 ) i d = l ，2 ，3 p ，o f 假设弹性体边界光滑：则c , j ( 尸) = 寺毛 利用边界元技术对积分方程( 2 2 4 ) 进行数值求解具体步骤如下： ( 1 ) 进行边界离散得r 的近似l ( 至少是线性的) ( 2 ) 对边界函数u , t 采用插值法，进行离散。 ( 3 ) 对边界节点按边界条件赋值，并建立节点信息。 ( 4 ) 以边界节点为配置点建立积分方程组。 ( 5 ) 把边界上的积分离散为边界单元上的积分求和，根据已知边界条件调整方 程中已知和未知量的位置，得出线性方程组。 ( 6 ) 求解线性方程组，得出全部边界量u ( 位移) 和t ( 面力) ( 7 ) 应用方程( 2 2 1 ) 求解原弹性问题( p 2 ) 。 本文采用四边形常单元对边界离散，假设边界被分为n e 个单元，对每个边界 节点p i ，( 2 2 4 ) 变为： z u ，= 0 q ( 2 2 5 ) j = ij = 1 其中：h f = j 玩( 只，o ，) a s ( q ，) + c 峨。l , m 2 l ，2 ，3 ( 2 2 6 ) g f = j u l m ( 只，q j ) d s ( q ，) ( 2 2 7 ) c = 1 2 l 2 当l j 时，蛾、g u 都没有奇异性，可用数值积分公式，本文采用九点g a u s s 数值 1 0 重庆大学硕士学位论文 2 三维弹性接触问题边界元解法 积分公式进行计算。 当捌时，h u 、g u 都具有奇异性，因此不能直接使用数值积分，对于g o ，由于使 用的是平面边界元，故可采用解析法，具体方法如下： g 扩= 肛( p ，9 ) d s ( 9 ) = 艋而1 。妒14 v ”言妾艄) ( 2 2 - 8 ) ， 徽分簟觉 “) 或 鬟2 且l 局蕊坐碌幕 ( h ) 图2 2 i ( a ) 积分单死( b ) 局部坐标糸 f i g2 2 1 ( a ) i n t e g r a le l e m e n t ( b ) l o c a lc o o r d i n a t e 设积分单元为四边形4 4 4 以( 如图2 2 1 ( a ) ) ，p 为形心，分别连接p a 。， p a ，p a ，p a 4 ，将其分成四个三角形a p a l a 2 ( 1 ) ，色蹦2 a ，( 2 ) ，甜驾a 4 ( 3 ) ， 朋。a ( a 4 ) ，规定每个三角形方向( 以右手螺旋为正) ，则： q t m = 喜占云南k 3 - 4 。) 瓦+ v m 谛( 9 ) ( 2 2 ” 以l 为例进行计算： 铲z ，丽1 7 1 【( 3 - 4 v ) 5 q + 毒- j l a s ( ”( 2 - 2 - 1 0 ) 作p h 上a j a 2 ( 见图2 2 1 ( b ) ) ，= p q ，令p = j 跗j ，记 b = ( 蹦，p q ) ， 0 2 = ( p h ，p q ) ，目= ( p h ，p q ) 。 在以点p 为极点，一p h 为极轴的极坐标系下，积分以的计算公式分别为： 扣斋惫f 。姗+ 丽两1e 严弓去，2 删目眩z ， 式( 2 2 1 1 ) 右端第二项中被积函数- o r 利用新的局部坐标系可以表达成在以点p 为 ( z 玩 极点，- f g 为极轴的极坐标系下的关于口的函数。 具体方法如下：将原坐标系旋转得到新笛卡尔坐标系。呱x ：x ：( 见图2 , 2 1 ( b ) ) ， 使得面平行于雨，瓦平行于而：，此时在新坐标系下有： 里堡垡羔翌型苎堑望鱼i 兰一 ! 三丝堂丝堡墅囹里垫墨重堡鲨 万o r = c o s o , 寿= s i n 口，寿= 。 眨z 国?弼础 假设矩阵爿一= ( ) 未为两个笛卡儿坐标系的变换矩阵，根据向量场梯度不变性有： 西 缸j 西 缸2 西 屯= i , , i 陲i 缸j f 口l l口1 2 = l 口：。口： 【吒，： 通过( 2 2 1 3 ) 式，_ o r ，_ o r ，- d r 表示成以口为变量的函数表达式，这样就可以 c x lo x ，c 在极坐标系下直接进行积分计算。 令：卜e 庐，跏口= p i n s e e 0 + 留嘶 ( 2 2 1 4 ) ，2 = 声；去) 2 ，跏p = 声( 口l l c 。s 口+ 口l ：s t n 咖d 目 2 【p s i n 口一2 口- - q ：p c o s 一十口磊p 1 n t g o + s e c o l - d 磊尸s i n 目1 2 ( 2 2 1 5 ) j i i2 而3 丽- 4 丽v ，1 + 而否而1l( 2 ，2 1 6 ) 同理 屯3 赤 慨1 _ o _ l r _ o r 嚣( 1 ) 2 而弓苦习乎棚庐( c 。s 口+ 口l ：s i n 颤0 2 ，c o s 口+ a 2 2s i n 目) 办 = 而而1 k u a ”p s i n o - ( a 1 1 c 2 2 + n 1 2 t 2 2 1 ) p c o s o + a ”a 2 2 p l n l t 9 0 + s e c o 一口。：口：p s i n 臼 l ： ( 2 2 1 7 ) j 1 3 、j 2 2 、j 2 ，、j 3 3 的计算依此类推，分别可以利用此方法得到解析积分表达式。 至于坐标系变换矩阵a “的计算则利用各节点在原坐标系下的坐标使用向量工 具即可直接计算。 计算过程如下： 仍然在朋。a ：中考虑，假设点p h ，和4 ：在原坐标系下的坐标分别为 ( x o ，y o ，z o ) ，( x i ，y i ，：1 ) ，( x 2 ，y 2 ，z 2 ) 。贝u 有： p h = ( 薯一x o ) i + ( y l j 0 ) j + ( ：。z o ) k h a 2 = ( x 2 一t v + ( n y i ) ，+ ( z 2 一z i ) 七 令： 1 2 )b 2 趴州， 试o c s v o o o = 重庆大学硕士学位论文2 三维弹性接触问题边界元解法 d l = j 、( x l x o ) 2 + ( m - y o ) 2 + ( 2 1 - z o ) 2 d 2 = 、( x 2 一x 1 ) 2 + ( 儿- y 1 ) 2 + ( z 2 一= 1 ) 2 则有： a = ( 一x o ) d l ；d 2 l = ( y 一y o ) d t ；a 叫= ( 毛一三o ) d l a 1 2 ；( 工2 一x i ) a 2 ；口2 z = ( y 2 一y 1 ) d 2 ；d 3 2 = ( z 2 一z 1 ) d 2 由于墓 = o ，因此无须计算，0 2 3 , 口3 ，。 对于h i 的计算，采用刚性位移的方法间接计算获得，其计算公式为： h “= - e 峨i = 1 , 2 n e ( 2 2 1 8 ) 篇 对n e 个节点循环计算后得如下方程组 h u = g t( 2 2 1 9 ) 其中h = ( h “) 3 n e x 3 n e g = ( g u ) 3 n e x 3 n e 根据边界条件将( 2 2 1 9 ) 中的已知量与未知量重新排列后可以写成 a x = f( 2 2 2 0 ) x 包含未知量( 面力，位移) 的列向量 2 3 接触问题的边界元法 考虑问题( p ，) ，对每个接触体分别可得边界积分方程( 2 2 4 ) ，一般地，可用 下式表示： c 矿( p ) “；( p ) + f 礞( p ，q ) “；( q ) d ( s 2 ) = f ( p ，q ) r ；( q ) d ( s ) ( 2 3 1 ) s s t k = a b i ，= 1 , 2 ，3 c f ( p ) = 寺岛 p ，q s 方程( 2 3 1 ) 任意单个积分方程利用2 2 节介绍的方法所得到的代数方程组都是 非线性的，因为接触区域边界上的位移和面力都是未知量。但是，方程( 2 t 3 1 ) 可 以通过弹性体之间在接触区域的面力与位移的相互关系( 即面力和位移接触条件) 来耦和以求达到减少整个问题未知量的目的。 2 3 1 接触条件 为了利用位移和面力接触条件，对接触区域作如下处理。 ( 1 ) 将两弹性体可能发生接触的边界区域以相同的方式划分。即具有相同的单 重庆大学硕士学位论文2 三维弹性接触问题边界元解法 元数，每个单元具有相同的形状。此时，来自两个弹性接触体的单元构成了单元配 对，节点则构成了节点配对。 ( 2 ) 引入局部坐标系( 见图2 3 1